Deciphering transcriptional regulators of banana fruit ripening by regulatory network analysis.

Fruit ripening is a critical phase in the production and marketing of fruits. Previous studies have indicated that fruit ripening is a highly coordinated process, mainly regulated at the transcriptional level, in which transcription factors play essential roles. Thus, identifying key transcription factors regulating fruit ripening as well as their associated regulatory networks promises to contribute to a better understanding of fruit ripening. In this study, temporal gene expression analyses were performed to investigate banana fruit ripening with the aim to discern the global architecture of gene regulatory networks underlying fruit ripening. Eight time points were profiled covering dynamic changes of phenotypes, the associated physiology and levels of known ripening marker genes. Combining results from a weighted gene co-expression network analysis (WGCNA) as well as cis-motif analysis and supported by EMSA, Y1H, tobacco-, banana-transactivation experimental results, the regulatory network of banana fruit ripening was constructed, from which 25 transcription factors were identified as prime candidates to regulate the ripening process by modulating different ripening-related pathways. Our study presents the first global view of the gene regulatory network involved in banana fruit ripening, which may provide the basis for a targeted manipulation of fruit ripening to attain higher banana and loss-reduced banana commercialization.

Introduction

Even though the growing world population, which is expected to reach 9.7 billion by the year 2050 (Gustavsson et al., 2011), renders providing for sufficient food supply a formidable task, a large fraction of horticultural crops is currently wasted, with losses occurring at many steps post‐harvest. Therefore, minimizing losses of horticultural crops and their products during the period from production to consumption is a central goal to meet current and future food demands (Sharma and Pongener, 2010; Suthar et al., 2019). Fresh fruits and vegetables are living tissues and are subject to continuous change and deterioration during the post‐harvest operations. Approximately 20%–30% of fresh horticultural products is lost during transport processes. This massive loss is caused by insufficient knowledge of proper handling practices, inadequate technology and/or poor storage equipment (Kumar and Kalita, 2017).

Fruit ripening is a highly complex biological process, and hence, precisely controlling it in order to avoid fruit spoiling would be greatly facilitated by a detailed understanding of the process. The ripening of fleshy fruit is characterized by significant changes of colour, texture, flavour, aroma and health‐promoting compounds, which collectively constitute the commercial quality traits relevant to consumers (Tucker et al., 2017). Previous studies have implicated that fruit ripening is controlled by the phytohormone ethylene and the expression of a large number of proteins and/or genes (Li et al., 2019; Ni et al., 2019). Of those ripening‐associated proteins, transcription factors (TFs) are of great interest as they play central regulatory roles in controlling the expression of various ripening‐associated genes and, thus, are responsible for phenotypic changes (Wang et al., 2019). For example, in tomato, TFs such as ripening inhibitor (rin; Vrebalov et al., 2002), colourless nonripening (cnr; Manning et al., 2006) and nonripening (nor; Giovannoni, 2004), together with their associated regulatory networks of influenced genes, have been identified as important cues in controlling fruit ripening. However, the regulation of fruit ripening is highly complex and intricate, in which multiple TFs and their corresponding target genes are involved (Lü et al., 2018). Thus, comprehensive and large‐scale studies of TFs and their interactions with their gene targets during fruit ripening are essential for developing an understanding of the molecular basis of fruit ripening. Ultimately, detailed molecular knowledge of the fruit ripening process may foster the development of novel biotechnological strategies for improving post‐harvest shelf life and fruit quality.

Experimental methods such as chromatin immunoprecipitation sequencing (ChIP‐Seq), yeast one‐hybrid screening (Y1H) and DNase hypersensitivity (DHS) assay have been successfully applied to study specific TF‐gene regulations. However, these methods are expensive, time‐consuming or technically demanding. Therefore, all such methods are of limited practicality and utility when implemented at a massive‐scale or genome‐wide analysis (Xiong et al., 2017).

With the rapid development of next‐generation sequencing technologies, transcriptome profiling via RNA‐Seq is a high‐throughput approach for generating a transcriptional map at the whole‐genome scale. Changes in transcriptomic data sets at multiple time points capture the progression of transcriptional changes through time and reveal different regulatory layers and pathways (Harkey et al., 2018). For example, using the RNA‐Seq technology, the regulation of maize leaf development has been well studied (Li et al., 2010; Wang et al., 2014; Yu et al., 2015), which provides useful information for developing C3 crops to equip with C4 photosynthesis modules, thereby increasing yield. As additional examples, the regulatory network of chickpea flower development and the dynamics of co‐expressed gene clusters in the gene transcriptional networks in Arabidopsis root have been revealed using time‐series transcriptomic data sets (Singh et al., 2013; Stringlis et al., 2018). These studies demonstrate that transcriptomic profiling is a powerful approach to discern regulatory networks and to discover new regulators involved in plant growth and development as well as response to abiotic and biotic stresses.

Bananas are the developing world’s fourth food crop in terms of gross value of production, after rice, wheat and maize. In 2016, about 17 million tons of bananas were globally traded, due to strong demand and growing health awareness of consumers worldwide (FAO, 2017). Banana fruits are a good source of minerals, vitamins, antioxidants, low glycemic carbohydrates and dietary fibre, which possess beneficial properties for human health. As a typical climacteric fruit, bananas are usually harvested and transported at pre‐climacteric, unripe stage. An artificial ripening is initiated by ethylene or ethylene‐releasing agents under controlled environments before marketing. However, the actual processes involved in ethylene regulating banana fruit ripening are not well understood, which makes a precise regulation of banana quality formation difficult. With the completion of banana genome sequencing (D’Hont et al., 2012), it has become possible to unravel the transcriptional regulatory networks of banana fruit ripening at the whole‐genome scale.

In this study, the pulp tissue of banana fruits treated with ethylene throughout the ripening process was sampled at eight time points. A Gene Expression Atlas was generated from RNA‐Seq data of banana pulp during the ripening process. A TF‐gene regulatory network was constructed based on WGCNA (weighted gene co‐correlation network analysis) and cis‐motif analysis. Based on the regulatory network, we identified critical TFs and cis‐elements involved in ethylene‐mediated fruit ripening of banana. Our findings provide valuable insights into the signalling and transcriptional processes involved in banana fruit ripening.

Results

Ethylene treatment accelerates banana ripening and induces ripening‐associated genes

Harvested banana fruits were treated with exogenous ethylene with untreated fruits serving as controls. Compared to control conditions, banana fruits exposed to ethylene treatment ripened more rapidly as assessed by phenotype, physiology evaluation and expression of ripening‐associated genes (Figure 1). Ethylene treatment resulted in prompting ripening phenotypic changes, whereas the control fruits showed almost no visual change (Figure 1a). Significant differences were also found for peel colour and firmness of bananas treated with ethylene. As shown in Figure 1b, the hue angle (h°) that represents the actual colour as perceived by an observer decreased markedly in peels of the treated fruit from 118° at Day 0 to 90° at Day 7 (Figure 1b), which indicates the peel colour varies from green to yellow. Similarly, the fruit firmness in ethylene‐treated bananas declined rapidly, with ~73 N at Day 0 decreasing to ~5 N at Day 7 (Figure 1b). Compared to the ethylene‐treated bananas, the peel colour and firmness of the control fruit remained unchanged. The ethylene‐treated bananas entered the climacteric phase at Day 3, showing the maximum level of ethylene production up to 3.62 µL/kg/h, while the control fruit retained low level of ethylene production during the entire storage period (Figure 1b). Further, we detected the expression of several known marker genes related to ethylene biosynthesis, cell wall modification, starch degradation and aroma formation in the process of banana ripening, which include MaACS1 (Liu et al., 1999), MaAMY3, MaSEX4 (Fan et al., 2018; Xiao et al., 2018), MaEXP1, MaEXP2 (Trivedi and Nath, 2004), MaOMT1 and MaMT1 (Guo et al., 2018). We found that expression of these ripening‐associated genes was significantly increased in the ripening stage of bananas exposed to exogenous ethylene, which is consistent with the exogenous ethylene production and texture softening. By contrast, no significant change in their expression was observed for the control fruit during the storage period tested (Figure 1c). These observations suggest that without ethylene treatment, banana fruits display no apparent change during the experimental period with regard to visual appearance, physiology and expression of ripening‐associated genes, but upon ethylene treatment undergo a series of physiological and molecular changes.

Figure 1

The variation of phenotype, physiology and expression of ripening marker genes of bananas with and without ethylene treatment in a 7‐day period post‐harvest. (a) Phenotypic variation of the natural ripening banana without ethylene treatment and ethylene‐induced bananas. (b) Variation of peel colour, firmness and ethylene production. (c) Variation of ripening‐association marker genes.

Gene expression changes of banana after ethylene treatment

In parallel with the performed ripening parameter assessments, a transcriptome profiling experiment using RNA‐Seq was performed on bananas at the same time points after ethylene treatment to capture the temporal gene expression changes during the ripening process. Three biological replicates were adopted in each time point. A total of 24 RNA‐Seq libraries for these samples were constructed and sequenced. After removing rRNA and low‐quality reads, the total length of clean reads ranged from 42 992 740 to 60 410 910 bp among the libraries, and 82.9% to 88.9% of the sequenced reads from each library could be perfectly mapped to the banana genome (Table S1). Hence, we conclude that the RNA‐Seq experiments have produced high‐quality transcriptomic data.

To verify the reproducibility and reliability of the RNA‐Seq transcriptome data, the expression of sixteen randomly selected genes was further investigated by qRT‐PCR analysis with gene‐specific primers (Table S2) in all samples of bananas treated by ethylene. As shown in Fig. S1, all sixteen genes showed a concordant expression pattern between the RNA‐Seq data and the qRT‐PCR results. Pearson correlation coefficients of RNA‐Seq and qRT‐PCR expression values ranged from 0.93 to 1.00 for all genes, supporting the reproducibility and reliability of the transcriptomic profiling data. Hence, the generated RNA‐Seq data set representing the ripening process‐associated gene expression variation proved to be highly reliable.

Principal component analysis (PCA) was carried out to analyse the relationships of the biological repeats in the samples as well as differences among samples. A total of 16 185 differentially expressed genes (DEGs) were subject to PCA (raw expression data referred to Data set S1). The DEG set was defined as all genes that were detected differentially expressed in at least one of the comparisons between ethylene treatment and control. The percentage of DEGs occupies 45.9% of the whole banana genome genes (the banana reference genome of Musa acuminate contains 35 276 genes). The result of PCA showed that three biological replicates of each time point clustered closely, indicating an acceptable variation within replicates at each time point (Figure 2a). Moreover, the PCA score plot also revealed a clear discrimination among the time points. As shown in Figure 2b, the samples could be separated into three groups. Group ‘one’ only contained the sample of 0 h (control), which clustered far from the remaining tested samples. This suggests that there are distinct patterns of gene expression between control and the ones after ethylene treatment. The transcripts have undergone an apparent change in the banana fruits after ethylene treatment. Group ‘two’ includes the samples clustered closely of 6 h, 12 h and 1 day. While Group ‘three’ were consisted of the samples of 3 day, 5 day and 7 day (Figure 2a). The differences of biological repeats and samples in Group ‘two’ and Group ‘three’ were quite small, indicating few differences in the time points of Group ‘two’ and Group ‘three’, respectively. A substantial transcriptional change was observed between Group ‘two’ and Group ‘three’, which meant that there was an apparent distinction 2 h to Day 1 and Day 3 to Day 7. The same trend was also observed in the phenotypic, physiological changes and gene expression evaluation of known marker genes (Figure 1).

Figure 2

PCA, number of DEGs and upset plot among different treatments. (a) PCA of three biological repeats of each time point after bananas treated with ethylene. (b) Number of DEGs of each ethylene treatment compared to control. (c) Upset plot (an alternative Venn diagram) of up‐regulated DEG of each comparison. (d) Upset plot of down‐regulated DEG of each comparison.

As illustrated in Figure 2b, the number of DEGs gradually increased following ripening, indicating an active dynamic machinery underlying fruit ripening at the transcriptional level. In addition, the number of down‐regulated genes is larger than that of up‐regulated ones, implying that genes with down‐regulation might play more important roles in banana fruit ripening than up‐regulated genes. Figure 2c and d indicates the overlapping genes and specific ones of up‐ and down‐regulated genes across the different time points. Day 3, Day 5 and Day 7 shared the largest number of overlapping genes, suggesting that they represent a stable state.

To further investigate the ripening process beyond general gene expression differences, transcription factors (TFs), which are likely critically involved in regulating fruit ripening, were specially investigated.

Expression patterns of TFs during the ripening process

A total of 1,315 TFs with different expression abundances were identified in the ripening process, which were classified into 55 families according to the PlantTFDB (Figure 3a, Jin et al., 2013). Figure 3 contains three types of information: the different expression TF number, the significantly over‐represented TF families in a certain period and the proportion of FPKM summary of a TF family relative to all TF genes. We calculated the number of differentially expressed TFs in each TF family during the ripening process. The significant analysis is determined by both whole gene number and the number of different expressional genes of a TF family, and the colour is determined by the expression of members of a TF family. A TF family could be white with its low‐expressed members, but it could be significant with lots of different expressional genes in its TF family members. Notably, the bHLH, ERF and MYB families of TFs comprised the largest numbers of differentially expressed TFs, with 108, 92 and 92 members, respectively. A significance analysis of the dynamic expression of related TF families across banana ripening stages was further generated. As observed before, at Day 3, the change in expression was most pronounced for TFs as well, supporting the notion that TFs have undergone a substantial transition at this stage of banana ripening. From the heatmap, some TF families have a higher proportion of FPKM than other TF families, such as bHLH, bZIP, C2H2, ERF, GARS, HD‐ZIP, MIKC‐MADS, NAC and TALE. In addition, the majority of TFs were down‐regulated (791 TFs), compared to 284 up‐regulated TFs (Figure 3b), which suggests that TF repressors might play more important roles in regulating fruit ripening. The temporal profiles ‘0’ and profile ‘1’ call for further attention since their variational trend fit the ripening process.

Figure 3

Expression dynamics of TF families. (a) The colour shows the total FPKM of all transcription factors (TFs) of a particular TF family. The asterisks indicate significantly over‐represented TF families at a time point. The first number in parentheses is the number of DEGs in a TF family, and the second number refers to the total number of genes in that TF family as annotated in PlantTFDB. (b) Expression profiles of all TFs with STEM (short time‐series expression miner). The log2(FPKM of ethylene‐treated sample + 1/FPKM of CK + 1) was used as input. Only 1200 TFs are shown as the genes were cancelled with fold change < 2 of the maximum input divided by minimum input per gene.

Although we have identified the TF families, which might play an important role during the ripening process, which members contribute to the ripening process exactly needs to be determined. Therefore, we further analysed the TFs as well as their target genes that were involved in banana fruit ripening.

Characterization of the critical TFs revealed by WGCNA and cis‐motif enrichment analysis

In order to identify the TFs, which potentially regulate the ripening process, the ripening‐associated genes should be characterized first. During the ripening process, the key physiological changes involve the ethylene production upsurge, starch degradation, aroma accumulation, texture softening and so on (Seymour et al., 1993). Therefore, the genes known to be involved in the processes of ethylene pathway (Jourda et al., 2014), starch degradation (Jourda et al., 2016; Xiao et al., 2018), sugar transport (Xiao et al., 2018), aroma biosynthesis (Asif et al., 2014) and cell wall modification (Asif et al., 2014) were selected as candidate ripening‐associated genes. Given that plant hormones have an indispensable role in fruit ripening (Kumar et al., 2014), genes encoding pathways of the following four phytohormones such as ABA, IAA, GA and cytokinin were also chosen. Taken together, those functional genes were made up of genes belonging to nine pathways: ethylene pathway, starch degradation, sugar transport, aroma biosynthesis, cell wall modification, abscisic acid (ABA) pathway, auxin (IAA) pathway, gibberellin (GA) pathway and cytokinin pathway.

The expression patterns of these functional genes are shown in Figure 4. It is apparent that the expression of many genes again varied dramatically from Day 1 to Day 3. Therefore, these ripening‐associated genes could be used as baits to identify key TFs regulating ripening process via correlation analysis.

Figure 4

Expression dynamics visualized as heat maps of nine ripening‐associated pathways and their associated genes. The left column to the right column indicates 0 h, 2 h, 6 h, 12 h, 1 day, 3 day, 5 day and 7 day after exposure to ethylene.

If a TF can regulate its gene targets directly, the TF‐gene targets probably have a high correlation relationship. Weighted gene co‐expression network analysis (WGCNA) can discern such relationships and generate an undirected weighted network with scale‐free topology. The network contains modules comprised by a cohort of genes with similar expression patterns in response to a given stimulus or physiological condition (Zaidi et al., 2019). The DEGs were subjected to a WGCNA. Results show that these DEGs are grouped into 6 modules, which are defined using colour codes. Each module contains genes with highly correlated expression behaviours (except the grey module with the left uncorrelated genes). Detailed information for each module is presented in Data set S2.

After the network was constructed by WGCNA, the connection of TFs and their target genes were preliminarily confirmed. Since the network was established on expression analysis alone, there could be many false‐positive connections. To try to reduce the rate of false‐positive connection, a motif enrichment analysis was carried out. If a TF binds and regulates a target gene directly, there should be the TFBS (transcription factor binding site) of this TF in the promoters of the target genes. TFBSs have been reported to be highly conserved in plants, indicating that their presence is of functional relevance (Franco‐Zorrilla et al., 2014; Lai et al., 2019). Therefore, a cis‐motif enrichment analysis was carried out. The results of WGCNA show that two modules, labelled as grey60 and lighbluesteel1 modules, were significantly associated with fruit ripening, because the two modules contained most of the ripening genes (324 of total 506 ripening genes, Table S3). The remaining modules were not further considered, since there were only few ripening‐associated genes in those modules. The 142 and 100 TFs, respectively, of the two selected modules, were further analysed with cis‐motif enrichment analysis. A total of 25 TFs along with their DNA‐binding motifs were identified, such as MaTALE1, MaMADS32, MaTALE2, MaERF54, MaMYB44, MaAS1, MaWRKY112, MaSPL1, MaZF‐HD8, MaERF83, MaC2H2‐4, MaG2‐1, MaZFP90, MaC2H2‐3, MaG2‐2, MabZIP100, MaC2H2‐5, MaMYB62, MaATAF, MaZFP48, MaNAC78, MabZIP6, MaZFP46, MaAP2‐46 and MabHLH8 (Figure 5). The ripening‐associated genes in the pathways regulated by the TFs are shown in Data set S3. Therefore, we predicted the hypothetical regulatory network of those TFs and ripening pathways. As seen in Figure 6, in addition to the TFs, which regulate the ripening‐associated genes, there were also many interconnections among the TFs. These findings suggest that fruit ripening is a coordinated process. Both, a TF and its interconnected TFs should be under consideration when regulating the fruit ripening process.

Figure 5

TFs targeting ripening‐associated pathway and associated potential binding site motifs.

Figure 6

Computed regulatory network of the 25 TFs and ripening‐associated pathways. The blue/red circles indicate down/up‐regulated TFs, respectively. Red edges highlight significant relationships among TFs. Blue edges denote relationships of TFs and their gene targets belonging to the nine ripening‐associated pathways.

Identification of cis‐motifs and gene regulation modules related to fruit ripening

To verify the reliability of the transcriptional regulatory network, we firstly utilized EMSA to test for the possible binding of TFs and their corresponding target genes. In total, 16 recombinant candidate TFs tagged with GST were successfully expressed and purified, which involves MaAP2‐46, MaATAF, MabZIP100, MaC2H2‐3, MaC2H2‐4, MaERF54, MaERF83, MaG2‐2, MaMYB44, MaMYB62, MaAS1, MaNAC78, MaSPL1, MaWRKY112, MaZF‐HD8 and MaZFP46. As expected, when the TF recombinant proteins were individually incubated with the biotin‐labelled probes derived from the promoters in corresponding targets, retarded mobility was clearly observed (Figure 7). These binding signals were out‐competed by the addition of excess unlabelled probes in the reactions, but not of the corresponding probe variants in which the core‐motifs were mutated (Figure 7), implying the specificity of the binding. Then, we performed yeast one‐hybrid (Y1H) assay to test direct binding of proteins to DNA sequences. The 16 segments containing the core‐motifs in promoters of target genes were cloned upstream of the Aureobasidin (AbA) resistance gene and transformed into Y1HGold strains as bait. After the careful selection of transformants on the SD/‐Ura plate and the examination of the minimal inhibitory concentration of AbA for the bait strains, the basal activities of 8 Y1HGold strains that were transformed with pAbAi‐MaPYL8pro, pAbAi‐MaPGMpro, pAbAi‐MaPEpro, pAbAi‐MaPL8pro, pAbAi‐MaXETpro, pAbAi‐MaEXPpro, pAbAi‐MaPL8pro and pAbAi‐MaXETpro were successfully inhibited, but they could grow well in the selection medium when transformed with effector plasmids expressing TFs (Fig. S2). On the contrary, the basal activities of the remaining 8 strains could not be inhibited by AbA, probably due to the technical limitation of this methods, and detection of their interaction was not suitable for Y1H assay. To further determine this association in planta, dual‐luciferase reporter assays in N. benthamiana leaves were performed. As shown in Figure 7, MaAP2‐46, MabZIP100 and MaC2H2‐3 activated the transcription of their target genes, while MaATAF, MaC2H2‐4, MaERF54, MaERF83, MaG2‐2, MaMYB44, MaMYB62, MaAS1, MaNAC78, MaSPL1, MaWRKY112, MaZF‐HD8 and MaZFP46 repressed their target genes’ transcription, suggesting the interaction between TFs and corresponding target genes. Additionally, due to the difficulties in banana transformation, we could not express the homologous genes in banana in a stable manner. It was found that a transient expression system in fruits could be used for investigating the functions of target genes and TFs (Wu et al., 2019; Xiao et al., 2018). Thus, we transiently overexpressed four TFs including MaZF‐HD8, MaC2H2‐4, MaATAF and MaAP2‐46 in banana fruit via Agrobacterium infiltration and determined the expression of their target genes. In accordance with the previous findings (Figure 7), the target genes of MaAP2‐46 were significantly elevated, while the target genes of MaZF‐HD8, MaC2H2‐4 and MaATAF were inhibited, when these TFs were successfully expressed in banana fruit (Fig. S3). Overall, these findings support the high reliability of the predicted transcriptional networks.

Figure 7

Validation of the binding of transcription factors to the cis‐elements in corresponding target genes by EMSA and tobacco‐transactivation assays. As for EMSA, GST‐tagged TF fusion proteins were incubated with cis‐element probes corresponding to different TFBSs and their mutated sequences (mTFBS) for EMSA. EMSA was conducted with labelled probes alone (lane a) or with combined purified TF and labelled probes (lane b). DNA‐binding specificity was tested by adding a 100‐fold excess of unlabelled probes (lane c) or a 100‐fold excess of unlabelled mutated probes (lane d). The core binding sequences are highlighted red. PWM, position weight matrix. For transactivation analysis, luciferase activities in tobacco leaves co‐infiltrated with the reporter and effector. The transactivation activity was determined by measuring the luciferase activities, with the activity in epidermal cells infiltrated with an empty effector construct defined as one.

Discussion

Construction of gene regulatory networks is a highly efficient approach to discern TF‐gene regulations in a variety of contexts, and WGCNA proved superior in identifying functionally associated gene pairs by examining pairwise gene expression (Xiong et al., 2017). By conducting transcriptomic analysis of pulp tissues from 0 to 7 days in ethylene‐induced fruit ripening, we established a co‐regulated gene network underlying banana ripening (Figure 6). Within the network, a total of 25 TFs together with their corresponding controlling pathways responsible for fruit ripening were identified. EMSA, Y1H and dual‐luciferase reporter system experiments indicated that 16 TFs, whose proteins could be successfully expressed and purified, were able to bind to their gene targets via the predicted cis‐motif elements (Figure 7), suggesting the validity of this network. Our work is the first report of uncovering the gene regulatory network associated with banana fruit ripening, which may have potential for improving storage life and quality traits of fleshy fruits through genetic engineering strategies.

It is generally accepted now that fruit ripening is regulated mainly by TFs. This was established, for example, based on three natural mutants in tomato such as rin, nor and cnr caused by reduced expression of MADS, NAC and SPL TFs, respectively (Giovannoni, 2004; Manning et al., 2006; Vrebalov et al., 2002). However, the number of TFs that were identified to control fruit ripening in a given plant is relatively small. Moreover, the relationships among these ripening‐regulated TFs have not been extensively explored. Therefore, it is urgent to develop an efficient method to unravel the regulatory relationships between the TFs and their gene targets, as well as the relationships among the TFs. To this end, a comprehensive analysis of dynamic transcriptome profiles is a powerful tool to decipher transcriptional regulatory relationships, because it reflects the correlated expression of pairs of TFs and their functional targets. Using this approach, we determined the gene regulatory network in banana fruit ripening. Among those TF nodes, most of them are interconnected to each other as well. Therefore, future studies in fruit ripening would aim at the modules of ripening‐regulated TFs.

A reliably predicted regulatory network should have a high significance for further gene function research. Therefore, we combined the published papers to evaluate the validity of the present research (Fasoli et al., 2018; Li et al., 2018). Among the TFs in the current regulatory model established here, some are homologous with senescence‐associated TFs identified in other species. For example, four TFs, namely MaMADS32, MaNAC78/MaATAF, and MaSPL1, which are homologues to tomato ripening regulators of RIN, NOR and CNR, respectively, were found in our network. As described previously, tomato RIN (Vrebalov et al., 2002), NOR (Giovannoni, 2004), and CNR (Manning et al., 2006) proved to be necessary for manifestation of fruit ripening. Other TFs, such as MaERF54/83, MaZF‐HD8, MaMYB44/62/MaAS1, MabHLH8, MaAP2‐46, MabZIP100, MaC2H2‐3/4/5 and MaWRKY112, are homologous genes of tomato SlERF6 (Lee et al., 2012), SlAP2a (Karlova et al., 2011), LeHB‐1 (Lin et al., 2008) and SlWRKYs (Wang et al., 2017), banana MaMYB3 (Fan et al., 2018), MabHLH6 (Xiao et al., 2018), MabZIP93 (Wu et al., 2019) and MaC2H2‐1/2 (Han et al., 2016a), which were reported to play roles in fruit ripening. Notably, some TFs including MaZFP46/48, MaTALE1/2 and MaG2‐1/2 were characterized to be involved in fruit ripening in this study for the first time. In general, ZFP TFs have been reported to be related to abiotic or biotic stresses (Mazumdar et al., 2017), while TALE TFs were shown to function in plant pathogen infection (Bogdanove et al., 2010), and G2‐like proteins constituted the GARP TF family with diverse functions of plant growth development and hormone signalling (Safi et al., 2017). The interrelationship among the TFs was complex. Further investigations focusing on the functions of these ripening‐associated TFs would greatly help decipher the molecular regulatory mechanisms of banana fruit ripening. The upstream TFs would be worth considered prior and verified experimentally before further usage.

TFs control gene expression through binding directly to cis‐regulatory elements, so‐called TF binding sites (TFBSs) in the promoters of their target genes. For this reason, the TFBSs provide valuable information to reveal the abundance and the temporal/spatial patterns of gene expression (Lu et al., 2019). In this context, using the MEME software, we identified 31 TFBSs in the regulatory model of banana fruit ripening (Figure 5). The data from our predictive TFBSs are in agreement with the available findings on binding specificity of corresponding TF families. Interestingly, transcription factors MaC2H2‐5, MaC2H2‐4, MaG2‐1 and MaATAF were found to bind to multiple different TFBS motifs, indicating that these TFs regulate the corresponding target genes during banana ripening in flexible manners. The characterization of the TFBSs in each TF would be fundamental to the understanding and construction of the model of transcriptional regulation in diverse plant biological processes.

We found that genes at 3 days after ethylene treatment displayed the most significant changes during fruit ripening (Figure 1). Combined with the ripening parameters tested, it seems plausible that formation of quality attributes of banana occurs mainly at 3 days after ethylene treatment. During this time point, we clearly observed the maximum level of endogenous ethylene production (Figure 1), which eventually contributes to significant changes in colour, aroma, flavour and texture in the pulp of banana. The climacteric rise of ethylene production is considered to be necessary to provide energy and intermediate substrates for biosynthesis of ripening‐associated compounds. Similarly, in tomato, the climacteric ethylene peak appears at the OR (orange) stage, when the quality characteristics (colour, size, shape, aroma, flavour and firmness) form (Liu et al., 2015). And we also further check the expression of the 25 TFs during the ripening period. As shown in Data set S4, all of the 25 TFs had an apparent expression modification on Day 3. It can be inferred that after ethylene treatment, the 25 TFs have undergone a fast adaptation and may orchestrate the gene expression changes that lead to the corresponding phenotype and physiology modifications that lead to fruit ripening.

In conclusion, we conducted WGCNA and cis‐motif analysis by using RNA‐Seq data from the pulp tissues of bananas treated by ethylene. A transcriptional regulatory model underlying banana fruit ripening was proposed, which consists of 25 TFs as well as 31 TFBSs in corresponding target genes responsible for hormone pathway, aroma formation, starch degradation, sugar transport and texture softening. Our findings not only offer a comprehensive overview of the banana fruit ripening network, but also narrow down the genes of interest in controlling banana ripening, which will provide valuable gene resource for molecular breeding strategy aiming at banana fruit quality improvement.

Methods

Plant materials and treatment

Mature green banana (Musa acuminate AAA group, Cavendish subgroup) bunches were sourced from farming fields in the Zengcheng district of Guangzhou, China (23°10′N, 113°61′E, alt. 97 masl.), and were subsequently transported to the laboratory within 2 h of harvest. Hands of banana were selected for uniformity of colour, size and weight, and absence of any physical injuries and diseases. Fruit hands were dipped into 0.05% SporGon (Decon Laboratories Inc., King of Prussia, PA, USA) for 3 min to inhibit pathogen infection, and dried at room temperature. They were randomly divided into two groups. The first group was treated by 100 µL/L ethylene for 18 h, while the control group was kept under constant regular air conditions for 18 h. After treatment, fruit hands were preserved in an environment with a temperature of 22°C and a relative humidity of 90% for ripening to take place. Samples were collected at 0 h, 2 h, 6 h, 12 h, 1 day, 3 day, 5 day and 7 day after ethylene treatment and immediately frozen in liquid N2 and then stored at −80°C for RNA extraction. The experiments were performed in three biological replicates.

Colour, firmness and ethylene production measurement

The ripening parameters in terms of colour, firmness and ethylene production were used to evaluate the ripening process of banana fruit. Peel colour was scored according to our previous report (Wu et al., 2019). Fruit firmness was determined as described previously (Han et al., 2016b). Exogenous produced ethylene content was measured according to Kuang et al. (2017).

RNA sequencing and data analysis

Total RNA was extracted from each pulp sample and stored at −80°C using the hot borate method (Wan and Wilkins, 1994). The library products were prepared for sequencing analysis via Illumina HiSeq™ 2500 (San Diego, CA). The raw reads were trimmed with Trimmomatic (Bolger et al., 2014) by discarding adapter sequences, reads containing more than 5% poly‐N and low‐quality sequences (reads with more than 50% percentage of low‐quality bases of quality value ≤ 5). Filtered reads were mapped to the banana reference genome Musa acuminate DH Pahang v2 as available at https://banana‐genome‐hub.southgreen.fr/download using Tophat v2.1.1 (Trapnel et al., 2012). Cuffdiff was applied to detect differentially expressed genes (DEGs) between control and ethylene treatment conditions. Significance values (P‐values) were corrected for multiple testing based on the false discovery rate (FDR) (Reiner et al., 2003). An FDR ≤ 0.05 and the absolute value of the log2 (fold change assessed by FPKM values) ≥1 were used as the threshold to determine significant differences in gene expression with FPKM meaning fragments per kilobase per million reads. FPKM of each treatment was derived from cuffdiff results. All raw reads were deposited to the National Center for Biotechnology Information Sequence Reads Archive (SRA) with accession number PRJNA598018.

Regulatory network construction by WGCNA (Weighted Gene Co‐correlation Network Analysis) and cis‐motif enrichment analysis

WGCNA was performed to generate the gene expression network based on correlation patterns among the DEGs. Genes were considered as DEGs only if their coefficient of variation between different treatments was less than 0.3 to eliminate noise. In total, 14 799 DEGs were used for analysis. The log2(FPKM of each treatment + 1) was used as input. Soft‐thresholding was set to 30 to make the network fit to a scale‐free topology. After the network was generated, the connection of transcription factors (TFs) and their regulatory target genes were preliminarily predicted. A sequence motif enrichment analysis was carried out for further to discard false‐positive connections. The network of module lightbluesteel was produced to cytoscape format with threshold 0.32, and the module grey60 was produced with threshold 0.20 to make the network suitable seize for further analysis. Based on the notion that if a TF could bind and regulate a target gene directly, the transcription factor binding site (TFBS) associated with this particular TF should be present in the promoter of the target genes. TFBSs have been observed to be highly conserved in plants (Franco‐Zorrilla et al., 2014; Lai et al., 2019), such that their detection can be regarded as informative. Motif detection in putative target genes was performed using MEME (Multiple EM for Motif Elicitation) with p‐value of 0.05, and TOMTOM was used to the compare the enriched motifs with known TFBS motifs to check if the enriched motifs belong to the TFBS reported for a particular TF with e value of 0.05. TFBS database information was obtained from JASPAR Core plants 2018 (http://jaspar.genereg.net), CIS‐BP (http://meme‐suite.org/db/motifs), DAP‐seq (affinity purification sequencing, O’Malley et al., 2016) and PBM (protein‐binding microarray, Franco‐Zorrilla et al., 2014). Only motifs known to occur in plants from those databases were used. The TF with connectivity more than 20 was analysed with cis‐motif analysis to avoid noises.

Quantitative real‐time reverse transcription polymerase chain reaction (qRT‐PCR)

The extracted RNA was used to produce cDNA by PrimeScript™ RT reagent kit with gDNA Eraser (TaKaRa, Shiga, Japan). Quantitative real‐time PCR (qRT‐PCR) was performed on a Bio‐Rad CFX96 Real‐Time PCR System (Bio‐Rad, Hercules, CA, USA) using the qPCR SYBR Green Master Mix (Promega, Madison, WI, USA). Gene expression levels were calculated according to the cycle threshold (Ct) value using MaRPS4 (ribosomal protein 4) as an internal control (Chen et al., 2011). Three replicate experiments were performed per sample.

Electrophoretic mobility shift assay (EMSA)

A total of 16 GST‐tagged TF recombinant fusion proteins were successfully expressed in Escherichia coli strain Rosetta (DE3) and then purified using the Protein Purification System (Clontech, Shiga, Japan). Fragments of length 60 bp containing the TFBS in the promoters of corresponding target genes were synthesized (Sangon Biotech, Shanghai, China) and biotin‐labelled using Pierce™ Biotin 3′ End DNA Labeling Kit (Thermo Fisher Scientific, Waltham, MA, USA). EMSA was carried using LightShift Chemiluminescent EMSA Kit (Thermo Fisher Scientific) according to the manufacturer’s directions. Different concentrations of non‐labelled actual and mutated probes were utilized for competition experiments.

Yeast One‐Hybrid (Y1H) assay

Y1H assays were performed using the Matchmaker Gold Yeast One‐Hybrid System kit (Clontech) as per the manufacturer’s guidelines. The fragments of target gene promoters were inserted into pAbAi vector to produce pAbAi‐bait plasmids which were subsequently linearized and transformed into yeast strain Y1HGold, and selected with a SD‐Ura plate. The ORF of TFs was ligated to pGADT7 to yield the construct AD‐TF. The constructs of AD‐TF were transformed into Y1HGold strain harbouring pAbAi‐bait and screened on a SD‐Ura/AbA plate. All transformation and screenings were performed twice.

Dual‐Luciferase Reporter (DLR) assay

DLR assays were carried out according to a protocol described by Xiao et al. (2018). 35S:REN served as a transformation efficiency control. LUC and REN activities were detected using Dual‐Luciferase Reporter Assay kit (Promega). The analysis was carried out using a Luminoskan Ascent Microplate Luminometer (Thermo Scientific) according to the manufacturer’s instructions. Six independent experiments (biological replicates) were performed.

Transient overexpression assay in banana fruit

Transient overexpression analysis in banana fruit was performed following the procedures described previously (Xiao et al., 2018). The ORF of candidate TFs was ligated to pCXUN vector under the control of a maize Ubiquitin promoter. Then, the recombinant constructs and the vector control were transiently expressed in banana fruit via Agrobacterium infiltration. Two or three days after infiltration, samples were collected for the measurement of gene expression. Experiments were performed in triplicate.

Conflicts of interest

The authors declare that they have no conflict of interests.

Author contributions

LC and WL conceived and designed the experiments. JK, CW and YG performed most of the experiments, and WS performed qRT‐PCR analysis. LC analysed the data. JK and LC wrote the manuscript. DW, JC and WL revised the manuscript. All the authors read and approved the manuscript.

Figure S1 Quantitative real‐time PCR validation of RNA sequencing results for ethylene‐induced ripening of banana fruit.

Figure S2 Validation of the binding of transcription factors to the cis‐elements in corresponding target genes by Y1H.

Figure S3 Transient overexpression of four candidate TFs in banana fruits and expression analysis of the target genes.

Table S1 Summary of transcriptome sequencing data.

Table S2 Primer sequences used in the present study.

Table S3 Modules deriving from WGCNA with TF numbers and ripening gene numbers.

Data set S1 The FPKM values for DEGs.

Data set S2 The modules of WGCNA.

Data set S3 The gene code of ripening pathway in Fig. 4.

Data set S4 The FKPM values of the 25 TFs in the network.

Click here for additional data file.