Genome-Wide Analysis of Basic Helix-Loop-Helix Transcription Factors to Elucidate Candidate Genes Related to Fruit Ripening and Stress in Banana (Musa acuminata L. AAA Group, cv. Cavendish)

The basic helix−loop−helix (bHLH) proteins are a superfamily of transcription factors (TFs) that can bind to specific DNA target sites, playing a central role in a wide range of metabolic, physiological, and developmental processes in higher organisms. However, no systemic analysis of bHLH TFs has been reported in banana, a typical climacteric fruit in tropical and subtropical regions. In our study, 259 MabHLH TF genes were identified in the genome of Musa acuminata (A genome), and phylogenetic analysis indicated that these MabHLHs could be classified into 23 subfamilies with the bHLHs from rice and Arabidopsis. The amino acid sequences of the bHLH domain in all MabHLH protein sequences were quite conserved, especially Arg-12, Arg-13, Leu-23, and Leu-79. Distribution mapping results showed that 258 MabHLHs were localized on the 11 chromosomes in the M. acuminata genome. The results indicated that 40.7% of gene duplication events were located in collinear fragments, and segmental duplications might have played a key role in the expansion of MabHLHs. Moreover, the expression profiles of MabHLHs in different fruit development and ripening stages and under various abiotic and biotic stresses were investigated using available RNA-sequencing data to obtain fruit development, ripening-specific, and stress-responsive candidate genes. Finally, a co-expression network of MabHLHs was constructed by weighted gene co-expression network analysis to elucidate the MabHLHs that might participate in important metabolic biosynthesis pathways in banana during development and the response to stress. A total of 259 MabHLHs were identified, and their sequence features, conserved domains, phylogenetic relationships, chromosomal distributions, gene duplications, expression profiles, and co-expression networks were investigated. This study systematically identified the MabHLHs in the M. acuminata genome at the genome-wide level, providing important candidate genes for further functional analysis. These findings improve our understanding of the molecular basis of developmental and stress tolerance in an important banana cultivar.


INTRODUCTION
Transcription factors (TFs) are important regulatory factors that interact with cis-elements to regulate the expression of genes in order to adapt to environmental stress in eukaryotes (Riechmann et al., 2000). Basic helix−loop−helix (bHLH) proteins contain a helix−loop−helix (HLH) structure (Kavas et al., 2016). The bHLH family is the second largest TF gene family (Feller et al., 2011) and plays an important regulatory role in the process of growth and development in plants (Pires and Dolan, 2010).
The typical bHLH domain usually contains 60 amino acids and can be divided into two functional regions: a basic area of 10-15 amino acids and a spiral area (HLH region) with 40 amino acids Feller et al., 2011). The alkaline region is located at the N-end of the BHLH domain, adjacent to the HLH region, and the region is mainly associated with the binding of TFs and specific sequences. The HLH region is located in the C-end of the bHLH domain, containing two adjacent α-helices. The homologous or heterogeneous dimers of the α-helices from the same or different bHLH TFs can be combined with motifs in the promoter of the target gene, thereby regulating the expression of the target gene .
An increasing number of bHLH superfamily members have been functionally characterized in various plant species since the first plant bHLH gene, which plays a key role in anthocyanin biosynthesis, was isolated from Zea mays (Ludwig et al., 1989). As an increasing number of genome sequences are being released, a variety of bHLH superfamily genes have been identified and analyzed in a wide range of plant species, such as Arabidopsis (Carretero-Paulet et al., 2010), soybean (Hudson and Hudson, 2015), wheat (Wang L. et al., 2019), Chinese jujube , rice (Li et al., 2006), peach , Chinese cabbage (Song et al., 2014), and cotton . However, the classification of a large number of TFs in the bHLH family has proved challenging (Carretero-Paulet et al., 2010). In animals, bHLH proteins have been organized into six groups, A through F, based on the conservation of specific amino acids and the bHLH and other domains (Atchley and Fitch, 1997). In plants, the classification of bHLH proteins is usually based on sequence homology within the conserved bHLH domain, and the number of subgroups ranges from 15 to 26 (Buck and Atchley, 2003;Pires and Dolan, 2010). Phylogenetic analyses of some atypical bHLH proteins have even extended that number to 32 subgroups (Carretero-Paulet et al., 2010).
In plants, most of the functionally characterized genes originate from Arabidopsis (Carretero-Paulet et al., 2010). The multifunctional bHLH TFs participate in a wide range of plant growth and developmental processes (Nesi et al., 2000;Ramsay et al., 2003;Groszmann et al., 2008;Kondou et al., 2008;Fujisawa et al., 2013) and secondary metabolism (Ludwig et al., 1989;Bailey et al., 2003;Heim et al., 2003;Zhang et al., 2012). bHLH TF family members also play key roles in the response to environmental stresses, such as drought (Seo et al., 2011), salt (Zhou et al., 2009), low temperature , and pathogen attack (Kim et al., 2012;Fan et al., 2014;Song et al., 2013). In banana, the study of bHLH is still in its infancy, and only a few bHLH genes have been functionally characterized (Peng et al., 2013;Zhao et al., 2013). The completion of the genome sequence of the Musa acuminata genome, a wild banana that belongs to subspecies M. acuminata ssp. malaccensis, has made a genome-wide analysis of bHLH genes possible (D'Hont et al., 2012). Banana (Musa. spp.) is an important staple food for many people in subtropical and tropical regions. It is rich in protein and carbohydrates, and banana production contributes significantly to many people's incomes (Aurore et al., 2009). Cavendish banana (M. acuminata AAA group cv. Cavendish) is a triploid (AAA) cultivar banana. However, the yield and quality of the banana fruits are often hindered by abiotic and biotic stresses, such as drought (van Asten et al., 2011), low temperature (Kang et al., 2007), salt (Xu et al., 2014), and several devastating diseases (Heslop-Harrison and Schwarzacher, 2007). One of the most important diseases, Panama disease or fusarium wilt of banana, is caused by the fungus Fusarium oxysporum f. sp. cubense (Foc) (Stover, 1962) and is widely regarded as one of the most destructive plant diseases in the world. Specifically, Foc tropical race 4 (Foc TR4) is a destructive disease of the commercial banana cultivar Cavendish (Ploetz, 2006;Visser et al., 2010). TR4 is responsible for most of the destruction to the commercial cultivar Cavendish, which is the main banana cultivar worldwide. Giant Cavendish Tissue Culture Variants (GCTCV) have acquired resistance to TR4 through somaclonal variation (Hwang and Ko, 2004). GCTCV-119 is the best Foc TR4-tolerant alternative cultivar for the Cavendish (Ploetz, 2015). Since the bHLH family has been proved to be associated with abiotic and biotic stress tolerance in many plant species, we systematically analyzed the bHLH gene family in banana to assess its potential relevance to fruit ripening and stress tolerance.

Identification MabHLH Genes in the Musa acuminata Genome and Phylogenetic Analyses
To analyze the genome-wide bHLH genes, we firstly downloaded all gene coding protein sequences of M. acuminata genome (A genome) DH Pahang v2 from the Banana Genome Hub 1 (Martin et al., 2016). Subsequently, We used the iTAK program 16 to identify TFs based on the consensus rules that are mainly summarized within PlnTFDB and PlantTFDB (Pérez-Rodríguez et al., 2009;Jin et al., 2013), and obtained all candidate MabHLH protein sequences. Finally, all candidate MabHLH protein sequences were further examined by the BLASTp and CDD 2 databases in NCBI. The bHLH protein sequences from rice and Arabidopsis were acquired from the RGAP 3 and TAIR 4 databases, respectively. The full-length MabHLHs protein sequence from M. acuminata, Arabidopsis, and rice were aligned using ClustalW.
Relationships were assessed using a bootstrap neighbor-joining evolutionary tree with 1,000 bootstrap replicates, created using MEGA 6.0 software (Tamura et al., 2013). The accession number of identified banana bHLHs is listed in Supplementary Table 1. The molecular weight and isoelectric points of the MabHLHs were predicted from the ExPASy database 5 . The sequence logo for bHLH domain was created by submitting the multiple alignment sequences to the WebLogo server 6 .

Chromosome Distribution and Gene Duplications
To determine the physical locations of bHLH genes, the starting and ending positions of all bHLH genes on each chromosome were obtained from the banana A genome database. MapInspect software was used to draw the images of the locations of the banana bHLH genes 7 .
Tandem and segmental duplications were also identified according to the plant genome duplication database (Lee et al., 2012). Examples of tandem duplication were identified based on physical chromosomal location: homologous MabHLH genes on a single chromosome, with no other intervening genes, located within 30 kbp of each other, were characterized as tandem duplication (Shiu and Bleecker, 2003;Du et al., 2013). Syntenic blocks were detected using MCSCAN (parameters: -a -e 1e-5 -s 5) , and all MabHLH genes located in the syntenic blocks were extracted. Circos (0.63) software was used to draw the images of the locations and synteny of the MabHLH genes 8 .

Plant Materials and Treatments
Cavendish banana (M. acuminata L. AAA group cv. Cavendish) is a triploid cultivar characterized by a high yield, high quality, and long-term storage. However, this variety is highly susceptible to Foc TR4. The banana fruits were obtained from the Banana Plantation of the Institute of Tropical Bioscience and Biotechnology (Chengmai,Hainan,20N,110E). For temporal expression analysis, developing banana fruits of 0, 20, and 80 DAF, representing fruit developmental stages of budding, cutting flower, and harvest stages, respectively, were collected from both banana varieties. BX fruits stored for 0, 8, and 14 DPH, representing the three progressive ripening stages based on color of the fruit, including green, yellowish green, and yellow, respectively, have been selected for post-harvest analysis. One month old banana seedlings were grown in Hoagland's solution (0.51 g/l KNO 3 , 0.82 g/l Ca(NO 3 ) 2 , 0.49 g/l MgSO 4 ·7H 2 O, 0.136 g/l KH 2 PO 4 , 0.6 mg/l FeSO 4 , 2.86 mg/l H 3 BO 3 , 1.81 mg/l MnCl 2 ·4H 2 O, 0.08 mg/L CuSO 4 ·5H 2 O, 0.22 mg/l ZnSO 4 ·7H 2 O, 0.09 mg/l H 2 MoO 4 ·4H 2 O) (pH 6.0) under greenhouse conditions (Arnon and Hoagland, 1939). The minimum and maximum temperatures in the greenhouse during the experiment were 25 • C and 30 • C, respectively, while the relative humidity oscillated between 55 and 80%. For salt and osmotic treatments, banana seedlings were irrigated with 300 mmol·L −1 NaCl and 200 mmol·L −1 mannitol for 7 days and havest the leaves without major vein to analysis. For cold treatment, banana plants were maintained at 4 • C for 22 h and havest the leaves without maine vein to analysis. For Foc TR4 treatments, we used the susceptible cultivar Cavendish banana and resistant cultivar GCTCV-119 (M. acuminata L. AAA group, cv. GCTCV-119) to inoculate Foc TR4 (Hwang and Ko, 2004). The roots of five-leaf stage banana seedlings were dipped in a Foc TR4 spore suspension of 1.5 × 10 6 condia/mL. The entire root system was then harvested at 0, 2, 4, and 6 days post-infection (DPI) (Wang Z. et al., 2012). All of the above samples were immediately frozen in liquid nitrogen and stored at −80 • C until RNA extraction for transcriptome analysis and QRT-PCR verification.

Transcriptomic Analysis
Total RNAs were isolated using a plant RNA extraction kit (TIANGEN, Beijing, China). Three µg of total RNA from each sample was converted to cDNA using a RevertAid First-Strand cDNA Synthesis Kit (Fermentas, Beijing, China). cDNA libraries were constructed using TruSeq RNA Library Preparation Kit v2, and were subsequently sequenced on an Illumina HiSeq 2,000 platform (San Diego, CA, United States) using the Illumina RNA-seq protocol. Two biological replicates were used for each sample. Gene expression levels were calculated as Reads Per Kilobases per Million reads (RPKM) (Mortazavi et al., 2008). Differentially expressed genes were identified with the read count of two replicates for each gene (fold change ≥ 2; FDR ≤ 0.001) (Audic and Claverie, 1997). A heat-map was constructed with MeV 4.9 and Java Treeview software according to the manufacturer's protocol. All the data of RNA-seq has been uploaded to deposit in the CNSA 9 of CNGBdb with accession number CNP0000292 (CNX0051086, CNX0051087, CNX0051090, CNX0051091, CNX0051097, CNX0051098, CNX0051099, CNX0051103, CNX0051104, CNX0051108, and CNX0051109) and analyzed as described in a previous study .

Weighted Gene Co-expression Network Analysis
Gene expression patterns for all identified genes were used to construct a co-expression network using WGCNA (version 1.47) (Langfelder and Horvath, 2008). Genes without expression detected in all tissues were removed prior to analyses. Soft thresholds were set based on the scale-free topology criterion employed by Zhang and Horvath (2005). An adjacency matrix was developed using the squared Euclidean distance values, and the topological overlap matrix was calculated for unsigned network detection using the Pearson method. Co-expression coefficients more than 0.55 between the target genes were then selected. Finally, we extracted the co-expression network of all MabHLHs and the network connections were visualized using cytoscape (Shannon et al., 2003).

Quantitative RT-PCR (QRT-PCR) Analysis
Total RNA was isolated from banana roots before treatment. First-strand cDNA was synthesized from 2 µg poly-(A) + RNA from each sample using AMV Reverse Transcriptase. SYBR Premix Ex Taq (TaKaRa, Dalian, China) was used in 25 µL reactions with 0.5 µL ROX reference dye. Primers (100 nM each) were mixed with the equivalent of 100 ng reversetranscribed RNA template per reaction. In all experiments, negative controls containing no template RNA were subjected to the same procedure to exclude or detect any possible contamination. Before proceeding with the actual experiments, a series of template dilutions were performed to determine the optimal template concentration necessary to obtain the maximal amplification of the target.
Each qRT-PCR was performed on a Stratagene Mx3000P (Stratagene, CA, United States) machine using SYBR chemistry. The thermal cycling conditions were as follows: 94 • C for 3 min, followed by 40 cycles of 94 • C for 15 s, 55 • C for 20 s, and 72 • C for 40 s. Reactions were performed in triplicate, and data were analyzed using MxProTM QPCR software (Stratagene, CA, United States). The MaActin transcript (Genebank accession numbers: EF672732) was used as a control. All the primers are listed in Supplementary Table 12, and the experiments were carried out with three biological replicates. The differences in C t values between the MabHLHs and MaActin transcripts were expressed as fold-changes related to MaActin.

Statistical Analysis
Statistical analysis was performed using Student's t-test. The experimental results obtained were expressed as the means ± standard deviation (SD). P values < 0.05 were considered statistically significant ( * ), and P values < 0.01 were considered highly statistically significant ( * * ).

Identification and Classification of MabHLHs in Banana
In total, 267 putative bHLH TFs were identified on the iTAK website (Zheng et al., 2016). The bHLH protein sequences encoded by non-representative transcripts were excluded. The remaining sequences were assessed for the existence of complete bHLH domains using the Conserved Domains Database (CDD). In total, 259 sequences were confirmed as putative members of the MabHLH family. The MabHLHs were randomly distributed on 11 chromosomes and were named MabHLH001 to MdbHLH258 based on their chromosomal locations. The highest number of MabHLH genes (29) was distributed on chromosome 3, while the least genes (14) were distributed on chromosome 11 (Figure 1)

Phylogenetic Analysis of the MabHLH Protein Family
To analyze the evolutionary relationships among the MabHLHs, 259 MabHLH proteins were aligned with 147 AtbHLH proteins from Arabidopsis (Li et al., 2006) and 157 OsbHLH proteins from rice (Pires and Dolan, 2010), and an unrooted phylogenetic tree was constructed by MEGA6 (Supplementary Table 2). The 259 MabHLHs were grouped into 23 subfamilies. In banana, there was no MabHLH member distributed in the subfamilies X, IIIX, and IVb. Furthermore, four bHLH proteins (MabHLH002, MabHLH003, MabHLH041, and MabHLH098) were not distributed in any subgroups, and were termed 'orphans' subgroups (Pires and Dolan, 2010). There were 45 MabHLHs distributed in subgroup XII, which was the largest; 25 MabHLHs distributed in subgroup Ia; and 20 MabHLHs distributed in the subgroup IIId + e and XV subfamilies. In these subgroups, the number of MabHLH members was far greater than that of Arabidopsis and rice, suggesting the significant expansion of MabHLHs in banana (Figure 2).

Conserved Amino Acid Residues in the MabHLHs Domain
Using CD-Search in NCBI and ClustalW align, all 259 protein sequences were found to contain the bHLH domain. The longest conserved domain was composed of 79 amino acids, including one alkaline region (basic) and one HLH region (helix 1) at the N-terminal, another HLH region (helix 2) at the C-terminal, and a loop located between helix 1 and helix 2, which was the least conserved (Figure 3). The frequencies of the 18 amino acid sites were higher than 50%. Four amino acid sites (Arg-12, Arg-13, Leu-23, and Leu-77) were conserved, and the frequencies of the conserved amino acids exceeded 95%. The frequencies of nine FIGURE 2 | Phylogenetic analysis of the bHLHs from Arabidopsis, rice, and M. acuminata genome. The neighbor-joining (NJ) tree was drawn using MEGA 6.0 with 1,000 bootstrap replicates.
amino acid sites (Glu-9, Arg-10, Arg-12, Arg-13, Leu-23, Pro-28, Ala-66, Leu-83, and Leu-79) were higher than 80%. All of the above amino acid sites were located in the basic, helix 1, and helix 2 regions. We found that the two helix regions were most conserved, followed by the basic region and then the loop region. Between the two helix regions, the most conserved site was located in helix 1 (Leu-23, 99%) (Figure 3).

Chromosomal Localization and Gene Duplication of MabHLHs
Genome chromosomal location analyses revealed that 258 MabHLHs were distributed on 11 chromosomes. As shown in Figure 1, the number of MabHLHs was irregular, although the MabHLHs were distributed on all of the 11 chromosomes. Although each of the 11 banana chromosomes contained MabHLHs, the distribution appeared to be uneven. There was no member located in the upper end of chromosome 2 (Figure 4). The distribution of MabHLHs on individual chromosomes also indicated certain physical regions with a relatively higher accumulation of gene clusters. For example, the MabHLHs on chromosomes 3, 6, 7, 8, 9, and 10 appeared to be congregated at the lower end and upper end of the arms (Figure 4).
Segmental duplication, tandem duplication, and retrotransposition are known to be key factors driving gene family expansion (Hurles, 2004;Freeling, 2009). Based on the FIGURE 3 | bHLH domain is highly conserved across all MabHLH proteins. The overall height of each stack represents the conservation of the sequence at that position. The black dot indicates the position of the 19 conserved amino acids previously identified by Atchley and Fitch (1997), and capital letters indicate over 50% conservation of amino acids among the 259 MabHLH domains. segmental fragment information in the M. acuminata genome (D'Hont et al., 2012), 92 colinear MabHLHs representing approximately 35.6% (92 of 258) of the total MabHLHs were found to be located in the syntenic blocks and had been segmentally duplicated (Figure 4 and Supplementary Table 3). Interestingly, we found that all the colinear MabHLHs within the syntenic regions belonged to the same group in the phylogenetic tree (Figure 2). By locating the genes on the chromosome, nine pairs were identified as tandemly duplicated (Figure 4). We did not find any retrotransposons. Based on our analysis, we propose that the expansion of MabHLHs was mainly via segmental duplication during the M. acuminata genome evolutionary process.

Expression Profile of MabHLHs During Fruit Development and the Postharvest Ripening Stage
To analyze the expression profiles of MabHLHs, RNA-Seq data were derived from the fruit development and ripening stages. We deleted 110 MabHLHs with Reads Per Kilobase of transcript per Million mapped reads (RPMK) values less than 1 in 0DAF, 20DAF, 80DAF_0DPH, 8DPH, and 14DPH. There were 97 and 95 MabHLHs highly expressed at 0 days after flowering (DAF) and 20 DAF (RPKM > 5), contributing about 37.5% and 36.7% of the total, respectively. We found that 51, 40, and 30 MabHLHs were highly expressed at 80 DAF-0 days post harvesting (DPH), 8 DPH, and 14 DPH (RPKM > 5), only contributing about 19.7, 15.4, and 11.6% of the total, respectively (

The Expression Level of MabHLHs in the Banana Seedlings in Response to Osmotic, Salt, Cold, and Foc TR4 Treatments
To analyze the expression profiles of the MabHLHs, RNA-Seq data were derived from banana seedlings in response to osmotic, salt, cold, and Foc TR4 treatments. To present the differentially expressed genes visually and exactly, we filtered out the genes with RPKM values less than 5 in both the control (0 day post infection) and treatments. Low-temperature treatment (4 • C) caused the 24 MabHLHs to be differentially expressed (log2 Ratio Cold/Control > 1), among which eight were up-regulated and 16 were down-regulated ( Figure 6A and Supplementary Table 6). Treatment with drought (200 mmol·L −1 mannitol) caused the 33 MabHLHs to be differentially expressed (log2 Ratio Osmotic/Control > 1), among which 12 were up-regulated and 21 were down-regulated ( Figure 6B and Supplementary Table 7). Treatment with salt (300 mmol·L −1 NaCl) caused 17 MabHLHs genes to be differentially expressed (log2 Ratio Salt/Control > 1), including five that were upregulated and 12 that were down-regulated ( Figure 6C and Supplementary Table 8). Infection with Foc TR4 resulted in the differential expression of 25 MabHLHs (log2 Ratio 2 DPI/0 DPI > 1), including 11 up-regulated and 14 downregulated genes (Figure 6D and Supplementary Table 9). These results suggest that these differentially expressed MabHLHs may play important roles in the response of banana to Foc TR4 infection. In the above osmotic, salt, cold, and Foc TR4 infection treatments, we detected 67 differentially expressed MabHLHs in banana, accounting for only 25.86% of the total. Generally, we found that most of the banana bHLH TFs were not involved in the response to various abiotic and biotic stresses. The 67 differentially expressed MabHLHs detected above may play important roles in the response to multiple stresses in banana.

Weighted Gene Co-expression Network of MabHLHs
Co-expression networks provide insights into the patterns of transcriptome organization and suggest common biological functions for networked genes. Weighted Gene Co-Expression Network Analysis (WGCNA) is a method frequently used to explore the complex relationships between genes and phenotypes. The distinct advantage is that WGCNA transforms gene expression data into co-expression modules, providing insight into signaling networks that may be responsible for phenotypic traits of interest (Horvath et al., 2006;Shi et al., 2010;Udyavar et al., 2013). To explore the functions, 259 MabHLHs were selected as "guide genes" to seek co-expressed genes using an RNA-Seq dataset from 11 different transcriptomes, including fruit development and ripening stages, banana seedling response to osmotic, salt, and cold treatment, and banana roots inoculated with Foc TR4 (Supplementary Table 4). A total of 4761 genes as the target node, whose expression patterns were closely correlated with 48 MabHLHs, were identified with weighted values larger than 0.5 (Supplementary Table 9). Following visualization using Circos, the co-expression network of MabHLHs was divided into three models. There were 30, 17, and three MabHLHs contained in models 1−3, respectively (Figure 7 and Supplementary Table 10).
According to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the 4761 genes, 18 KEGG pathways were enriched, mainly including biosynthesis of secondary metabolites (ko01110), plant hormone signal transduction (ko04075), plant−pathogen interaction (ko04626), and phenylpropanoid biosynthesis (ko00940) (Figure 8 and Supplementary Table 11). Notably, the differentially expressed genes in banana in response to Foc TR4 infection, including MabHLH032, MabHLH062, MabHLH066, MabHLH067, MabH LH079, MabHLH153, MabHLH178, MabHLH185, MabHLH186, MabHLH210, and MabHLH252, were included in the weighted co-expression network of the MabHLHs. These results demonstrate that these MabHLHs might be important in regulating the response of banana to Foc TR4 infection.

Expression Patterns of MabHLHs During the Interaction of Banana Plantlets With Foc TR4
We discovered that there were 11 MabHLHs involved in the above weighted co-expression network that were differentially expressed during the interaction of banana with Foc TR4, and all of them decreased significantly at 2 DPI in the RNA-Seq data. MabHLH062, MabHLH066, MabHLH067, MabHLH079, MabHLH153, and MabHLH185 were highly expressed (RPKM value > 20) at 0 DPI (Supplementary Table 9).
These six differentially expressed MabHLHs were selected for quantitative real-time (qRT)-PCR analysis of their expression patterns in the response of the resistant (GCTCV-119) and susceptible (Cavendish) banana cultivars to Foc TR4 infection. RNA was extracted from the roots of two cultivars at 2 DPI, 4 DPI, and 6 DPI and subjected to qRT-PCR analysis. A mock treatment (0 DPI) was carried out using Hoagland's solution as a control. In the Cavendish banana, all four MabHLHs decreased significantly (P < 0.05, the same below) at 2 DPI, 4 DPI, and 6 DPI. We also found that six MabHLHs showed the same trends and consistent results between the RNA-Seq data and qRT-PCR data at 0 DPI and 2 DPI (Figure 9). These results indicated that the RNA-Seq data were suitable for detecting the expression patterns of MabHLHs. Moreover, in the GCTCV-119 banana, MabHLH062, MabHLH067, MabHLH153, and MabHLH185 increased significantly at 2 DPI and reached the highest relative expression levels of 10.98-fold, 9.88-fold, 4.59fold, and 4.97-fold at 6 DPI, respectively. MabHLH066 was not differentially expressed at 2 DPI and 4 DPI, and only decreased markedly at 6 DPI. MabHLH079 was also not differentially expressed at all the time points (Figure 9). In summary, these results indicated that MabHLH062, MabHLH067, MabHLH153, and MabHLH185 were involved in banana resistance to Foc TR4 infection in GCTCV-119 cultivar.

DISCUSSION
Whole-genome identification of bHLH genes may provide global insight into the basis for the functional diversification of these genes during evolution. The bHLH TFs comprise a large family in higher plants, and numerous studies have shown that bHLH TFs are involved in diverse biological processes in plant growth, development, and stress responses . In plants, bHLH is the second largest family of TFs and includes 319 members in soybean (Hudson and Hudson, 2015), 230 in Chinese cabbage (Song et al., 2014), 167 in rice (Li et al., 2006), and 162 in Arabidopsis (Carretero-Paulet et al., 2010). In our study, we identified 259 MabHLHs from the M. acuminata genome, which has been the second-highest number of bHLHs detected, to date. The characteristics of the number of encoding amino acids, isoelectric point, and molecular weight of the MabHLHs were similar to previous results in Arabidopsis and other plants (Carretero-Paulet et al., 2010).
Previous studies have shown that over 50% of plant bHLH proteins contain five alkaline amino acid residues in the conserved domain, as well as a highly conserved HER motif (His5-Glu9-Arg13), of which the predicted function is to bind with the target DNA motif (Atchley and Fitch, 1997;Massari and Murre, 2000;Toledo-Ortiz et al., 2003). We also found that all MabHLH proteins have conserved HER motifs. Compared with the conserved domain sequences of the HLH region from FIGURE 8 | KEGG pathway enrichment analysis of the genes in the co-expression network of the MabHLH genes. A total of 4809 (4761 target node genes and 48 MabHLHs) genes were used in the enrichment analysis. The gene set enrichment was analyzed using hypergeometric testing. The Q-value was calculated using the FDR (false discovery rate) adjustment method for correcting multiple hypothesis testing. The top 20 pathways are shown. Q-values represent the significance of the enrichment. Circles indicate the target genes, and the size is proportional to the number of genes.
Arabidopsis (Feller et al., 2011), rice (Li et al., 2006), tomato (Sun et al., 2015), and Chinese cabbage (Song et al., 2014), we found a similar frequency of amino acid residues, and the amino acid residue sites in  were highly conserved in the HLH region of the MabHLH proteins (Figure 3). We also found that the conservation of Helix1 was higher than that of Helix 2, which corroborates previous results in Arabidopsis (Feller et al., 2011), rice (Li et al., 2006), and tomato (Sun et al., 2015). Therefore, these findings indicated that the bHLH domain in the MabHLH proteins was highly conserved, and all members of MabHLHs are typical bHLH genes from the M. acuminata genome.
The M. acuminata genome has undergone three episodes of genome-wide duplication (WGD) events (α-, β-, and γ-WGD), and the duplicated segments included in the Musa ancestral blocks cover 222 Mb and contain 26,829 genes FIGURE 9 | Expression patterns of six MabHLH genes in Cavendish and GCTCV-119 inoculated with Foc TR4 by qRT-PCR. The results are presented as differential relative transcript abundance. The data represent the mean ± standard deviation (SD) of three replicates. The y-axis shows the transcript fold-change relative to that in the control (0 DPI). * and **significantly different from the control (0 DPI) at P < 0.05 and 0.01, respectively. (D'Hont et al., 2012). We found that 92 MabHLHs were segmentally duplicated genes, and nine pairs of MabHLHs had been tandemly duplicated (Figure 4), which suggested that segmental duplication was the main pathway of MabHLH expansion in the M. acuminata genome.
Plant bHLHs play important roles in growth and development (Ding et al., 2009;Ikeda et al., 2013), stress tolerance (Kiribuchi et al., 2005;Miura et al., 2007;Long et al., 2010;Song et al., 2013;Fan et al., 2014), and signal transduction (Kiribuchi et al., 2005;Seo et al., 2011). In banana fruits, five MabHLHs were found to participate in methyl jasmonate (MeJA)-induced chilling tolerance (Peng et al., 2013), and MaMYC2 belonging to the bHLH gene family can bind with MaICE1 to regulate induced chilling tolerance in fruit (Zhao et al., 2013). In banana, MabHLH6 was increased both in transcript and protein levels, and positively regulated the expression of starch degradationrelated genes during fruit ripening process . MabHLH6 was negatively regulated by MaMYB3, which delayed the ripening of banana fruit (Fan et al., 2018). The transcript level of MabHLH7 was obviously increased during banana fruit ripening, which was ethylene-inducible and nuclear-localized (Song et al., 2020). In our results, the expression of MabHLH 024 (same with MabHLH7) and MabHLH139 (same with MabHLH6) was also increased during 80DAF_0DPH, 8DPH and 14DPH (Figure 5), which is consistent with MabHLH6 Song et al., 2020). Furthermore, about 40% of MabHLHs were highly expressed during banana fruit development (0 DAF and 20 DAF). The expression profile results indicated that the MabHLHs were significantly involved in fruit development. Banana fruit is a typical climacteric fruit that releases a large amount of ethylene during the postharvest ripening process (Liu et al., 1999). We found that most of the MabHLHs were lowly (RPKM < 5) or not expressed in the ripening stage, indicating that most of the MabHLHs may not be involved in the postharvest ripening of banana fruit. At the same time, we found that MabHLH229, MabHLH035, MabHLH152, MabHLH263, MabHLH155, and MabHLH157 were differentially expressed and gradually decreased at 0 DPH, 8 DPH, and 14 DPH, suggesting that these genes might play negative regulatory roles in banana fruit during the postharvest ripening stage.
Many studies have shown that bHLH is involved in the tolerance to abiotic stress. In rice, OsbHLH006 and OsbHLH148 participate in the traumatic response and drought stress through regulating the expression of genes related to the rice jasmonic acid signaling pathway (Kiribuchi et al., 2005;Seo et al., 2011). In Arabidopsis, AtICE1 TFs participate in plant stomatal development and low temperature stress (Miura et al., 2007;Kanaoka et al., 2008). In our results, we found that more than 20% of genes were differentially expressed in the cold, drought, and salt stress treatments. These results indicate that MabHLHs play important roles in the banana response to various abiotic stresses, such as drought, high salinity, and low temperature.
The bHLH TF family also plays an important role in plant disease resistance. In rice, OsRAI1 encodes a putative basic helix-loop-helix TF. The expression of the OsRAI1 gene can be induced by pathogens and can promote the expression of OsPAL1 and OsWRKY19 to participate in the defense response to rice blast (Kim et al., 2012). In Arabidopsis, AtHBI1 (homolog of BEE2 interacting with IBH1OMOLOG OF BEE2 INTERACTING WITH IBH1) can form a basic complex with the PRE gene and IBH1 gene (PRE-IBH1-HBI1), playing a regulatory role in plant resistance to pathogen infection and growth (Fan et al., 2014). Members of the third bHLH TF subgroup, which includes AtbHLH3, AtbHLH13, AtbHLH14, and AtbHLH17, act as transcription repressors to interact with jasmonate ZIM-domain proteins to repress jasmonate responses (Song et al., 2013). In our results, we found that six MabHLH genes were highly expressed (RPKM > 20) in the banana roots during Foc TR4 infection (Figure 9 and Supplementary Table 9). Traditional biological research starts with a single gene or protein, which can only explain part of the system and makes it difficult to comprehensively explore the whole system. Gene networks are valuable for gene function discovery and candidate gene prioritization. Compared with other methods, WGCNA can better preserve the characteristics of biological networks, and the modules can better reflect the relationship among functions and different biological processes (Kadarmideen and Watson-Haigh, 2012). In the co-expression network, the plant−pathogen interaction (ko04626) pathway was enriched, indicating that MabHLHs might be involved in the banana response to Foc TR4 infection. We also found that 11 differentially expressed MabHLHs were contained in the co-expression network. In our study, the relative expression values of MabHLH062, MabHLH067, MabHLH153, and MabHLH185 were larger in the resistant cultivar than in the susceptible cultivar (Figure 9). These results indicated that these MabHLHs were involved in banana resistance to Foc TR4 infection and the mutation site in GCTCV-119 cultivar selected by somatic variation might directly or indirectly affect the expression of these MabHLHs in the banana response to Foc TR4 infection.

CONCLUSION
We systematically identified 259 MabHLHs in the M. acuminata genome. A total of 258 MabHLHs were located on 11 different chromosomes and could be classified into 23 main groups based on phylogenetic analysis with bHLHs from Arabidopsis and rice.
We analyzed the expression profiles of the MabHLHs at different stages of fruit development and ripening and found numerous MabHLHs that participate in the banana response to abiotic and biotic stresses. Finally, the co-expression network of MabHLHs was constructed using WGCNA to elucidate the MabHLHs that might participate in important metabolic biosynthesis pathways in banana. The qRT-PCR results suggested that MabHLH062, MabHLH067, MabHLH153, and MabHLH185 were involved in banana resistance to Foc TR4 infection. This comprehensive study improves our understanding of the MabHLHs associated with fruit development, ripening processes, and stress responses and will establish a foundation for future studies on genetic improvement in banana.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the data has been uploaded to deposit in the CNSA (https://db.cngb.org/ cnsa/) of CNGBdb with accession number CNP0000292.

AUTHOR CONTRIBUTIONS
ZW, BX, and ZJ conceived the study. ZW, CJ, J-YW, J-HL, CC, and H-XY analyzed the bHLH gene's data from M. acuminata genome and performed the bioinformatics analysis. ZW, CJ, CC, and H-XY performed the Foc TR4 inoculation, RNA isolation, and qRT-PCR. ZW wrote the manuscript. All authors read and approved the final manuscript.