Genome-Wide Identification and Characterization of bHLH Gene Family in Hevea brasiliensis

1. Introduction

The rubber tree (Hevea brasiliensis) is the primary source of natural rubber, which is an indispensable raw material for industrial production. More than 90% of natural rubber in the world is produced from rubber trees[1]. The regulation of rubber biosynthesis and related physiological processes is governed by complex genetic networks in which transcription factors play a crucial role. Among these, the basic helix-loop-helix (bHLH) transcription factor family has recently gained attention for its involvement in various biological processes in plants, including growth, development, and responses to environmental stress.

The basic helix-loop-helix (bHLH) family is the second-largest family of eukaryotic transcription factors, following the MYB family, and is known for its basic helix-loop-helix (HLH) domain, which consists of approximately 60 amino acids, including a basic region and an HLH region. The basic region, located at the N-terminus, is rich in basic residues and serves as a DNA binding site. The HLH region, at the C-terminus, acts as a dimerization domain, enabling the formation of protein dimers through specific interactions, and is primarily composed of hydrophobic amino acids forming two amphipathic α-helices separated by a loop segment with a variable sequence and length[2].

bHLH transcription factors were found in both animals and plants. In animals, they are typically classified into six groups (A to F) based on phylogenetic relationships and DNA-binding patterns[3,4]. However, there are differences in the classification of plants and animals because of the relative independence of bHLH transcription factor genealogy[5]. To date, no definitive classification standard for bHLH transcription factors has been established in plants. According to bHLH protein sequence homology, conservation of specific amino acids at key positions, and the presence of other conserved functional domains outside the bHLH domain, researchers have classified bHLH transcription factors of plants into 15-35 subfamilies (Figure 1). For instance, in Arabidopsis, 162 bHLH proteins have been characterized and grouped into 12 major groups and 25 subgroups[6,7]. Lorenzo Carretero-Pauletd, by including lower plants such as algae and mosses in the analysis, has classified the bHLH transcription factors into 32 subgroups from a broader evolutionary perspective[8]. Rice (Oryza sativa L.) has 167 bHLH family members, which are classified into 22 subfamilies[9]. Tomato (Solanum lycopersicum) possesses 152 bHLH genes, organized into 26 subfamilies, whereas[10] cucumber (Cucumis sativus L.) has 142 bHLH genes distributed among 32 subfamilies[11]. In potato (Solanum tuberosum), 124 bHLH members have been recognized and divided into 15 subfamilies[12].

The ubiquitous presence of bHLH transcription factors further underscores their profound genetic significance and functional specialization throughout the plant kingdom. Recent studies have shown that the influence of the bHLH superfamily extends across a broad range of plant growth and developmental processes. For instance, bHLH transcription factors are instrumental in regulating a multitude of biological processes, such as seed germination[13,14], flowering time control[15,16,17], responses to iron deficiency[18,19], mechanisms for cadmium tolerance[20], and enhancement of cold tolerance[21]. Furthermore, they play a pivotal role in mediating responses to light signals[22,23,24] and phytohormones, notably jasmonic acid (JA) [25,26] and abscisic acid (ABA) [27,28,29], which are vital for environmental stress adaptation and physiological regulation. Collectively, these roles emphasize their critical importance in the intricate networks governing plant stress responses and development.

Prior to this study, researchers analysed the functions of various bHLH factors in H. brasiliensis. For instance, Yamaguchi et al. subjected four HbbHLH genes to RNA-seq and in vitro binding assays, discovering that HbMYC2 might function as a receptor in JA signalling, contributing to latex cell differentiation. In contrast, the other two HbbHLH genes likely promote rubber biosynthesis[30]. Chen et al. found that one bHLH transcription factor, ICE-like transcription factor HbICE2, is involved in the cold stress tolerance mechanism regulated by jasmonic acid in rubber trees[31]. Guo et al. found that HbMYC2b significantly increased the activity of the HbSRPP promoter, suggesting that HbMYC2b may be a positive regulator of HbSRPP in rubber trees[32]. However, our understanding of the overall function and characteristics of the bHLH family genes in rubber trees remains limited. In this study, we aimed to identify and characterize the complete bHLH family members of the rubber tree and investigate their transcriptional changes during different treatments, thereby providing clues for further functional characterization of genes.

2. Materials and Methods

2.1. Identification of Hevea bHLH Genes

The genome sequences of H. brasiliensis originated from a preliminary assembly (unpublished) performed by our research group, and the Hidden Markov Model (HMM) profile of bHLH (PF00011) was downloaded from the Pfam protein database (http://pfam-legacy.xfam.org/, accessed on March 25, 2024)[33]. Putative HbbHLH sequences were identified using the profile HMM, generated, and calibrated using the HMMER[34] software with default parameters. Then, the candidates were checked for the conserved bHLH domain using Batch CD-search (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi, accessed on 25 March 2024) [35], Pfam[33], and SMART (https://smart.embl.de/, accessed on 25 March 2024) [36]. The remaining sequences were preserved as the final collection of the putative HbbHLH sequences. The amino acid count, molecular weight, isoelectric point, instability index, aliphatic index, and grand average of hydropathicity (GRAVY) of the proteins encoded by the HbbHLH genes were determined using TBtools[37]. Subcellular localization predictions were made using the WoLF PSORT[38] web server (https://wolfpsort.hgc.jp/, accessed on 25 March 2024). All domain sequences of HbbHLH genes were aligned with ClustalW[39], and the alignment results were submitted to the online website WEBlogo (https://weblogo.berkeley.edu/logo.cgi, accessed on 25 March 2024)[40] to create a sequence logo.

2.2. Sequence Conservation and Phylogenetic Analysis

All domain sequences of HbbHLH were aligned with ClustalW[39], and the alignment results were submitted to the online website WEBlogo (https://weblogo.berkeley.edu/logo.cgi, accessed on 25 March 2024) [40] to create a sequence logo.

The bHLH amino acid sequences of A. thaliana were obtained from The Arabidopsis Information Resource (TAIR, https://www.Arabidopsis.org/, accessed on 26 March 2024). The AtbHLH and HbbHLH sequences were aligned using ClustalW[39] with default parameters, and the resulting alignments were subsequently used to construct an unrooted neighbor-joining phylogenetic tree using MEGA7[41], employing the following parameters: 1000 bootstrap replicates, Poisson model, and pairwise deletion of gaps. Thereafter, the phylogenetic tree was visually rendered using evolview (https://www.evolgenius.info/evolview/#/treeview, accessed on 27 March 2024) [42].

2.3. Gene Structures, Conserved Motifs, Promoter and Chromosomal Location Analysis

The structures of the HbbHLH genes were analysed using TBtools[37]. Conserved motif analysis was conducted with the MEME Suite[43] online program (https://meme-suite.org/meme/tools/meme, accessed on 30 March 2024), setting the maximum number of motifs to 20. Additionally, the 1500 bp upstream DNA sequences of the HbbHLH genes were submitted to the PlantCARE database (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/, accessed on 27 March 2024) [44] for the prediction of cis-acting regulatory elements, with the results visualized using TBtools[37]. HbbHLH genes were mapped to chromosomes using MapChart software[45].

2.4. Gene Duplication Patterns Analysis

The entire protein sequence of the H. brasiliensis genome served as queries in a self-alignment process executed with Blastall, employing an E-value threshold of 1e-10. The results were analysed using MCScanX[46] to produce segmental duplications and tandem duplications across the whole genome. Segmental duplications belonging to the bHLH family were extracted to draw a collinearity map using CIRCOS[47].

2.5. Gene Ontology Enrichment Analysis

To investigate the functional enrichment of the 180 HbbHLH genes, Gene Ontology (GO) analysis was performed based on gene function annotations obtained from previous studies within the research group. The R package ‘clusterProfiler’[48] was employed for this analysis, with significance thresholds set at Q < 0.05 and P < 0.05.

2.6. Interaction Network of Hevea Homologues in Arabidopsis

Homologous genes of H. brasiliensis were identified in A. thaliana using OrthoVenn3 (https://orthovenn3.bioinfotoolkits.net/, accessed on 30 August 2024)[49]. The HbbHLH homologous gene sets were further analyzed in STRING (https://cn.string-db.org/, accessed on 30 August 2024) with a minimum required interaction score of 0.4, limiting the results to no more than 50 interactors. Genes with homologues in Hevea were filtered from the interaction results and visualized using the Cytoscape[50] software.

2.7. Expression Patterns Analysis

The expression levels of the HbbHLH genes were quantified using the average of normalized expression data (measured in FPKM—Fragments Per Kilobase of transcript per Million mapped reads) obtained from HeveaDB (http://hevea.catas.cn, accessed on 22 May 2024)[51]. An expression heatmap was generated using the R package ‘pheatmap’ (https://CRAN.R-project.org/package=pheatmap). Expression pattern clustering was performed using Muffz[52].

The rubber tree cultivar RY7-33-97, cultivated at the Experimental Farm of the Chinese Academy of Tropical Agricultural Sciences (19°28′ N, 109°29′ E), was treated individually with 1% ethrel and methyl jasmonate. Five tissue samples (bark, leaves, latex, female flowers, and male flowers) were collected without prior treatment. Latex samples were collected at predefined post-treatment intervals: immediately preceding the treatment (0 h) and subsequently at 4, 8, 12, 24, and 48 h. For the branch treatments, four treatments were applied, and rubber tree leaves were sampled at 0, 2, 8, and 24 h post-treatment. Three biological replicates were used for each treatment.

Total RNA was extracted using an RNA Prep Pure Plant Plus Kit (TIANGEN, Beijing, China). The PrimeScript™ RT Reagent Kit with gDNA Eraser (Takara, Dalian, China) was used for first-strand cDNA synthesis. RT-qPCR was conducted using the2 × Q3 SYBR qPCR Master Mix (TOLOBIO, China) with primers designed using Primer3(Table S5), and the rubber tree YLS8 gene (Gene ID: 110644482) was used as an internal reference gene. The relative expression levels were calculated using the 2-ΔΔCT method. Each analysis was repeated in at least three independent experiments.

3. Results

3.1. Genome-Wide Identification of Hevea bHLH

Using Hidden Markov Model search results and screening through the SMART and CDD databases, 180 rubber tree bHLH family genes were identified and named HbbHLH1-HbbHLH180 based on their chromosomal locations. These genes display distinct features (Table S1); all HbbHLH genes have negative Grand Average of Hydropathicity values, indicating that they are hydrophilic. Subcellular localization analysis revealed that the vast majority (91%) of HbbHLH genes were predicted to localize to the nucleus, highlighting their potential roles in transcriptional regulation, while the remaining genes were localized to other organelles such as mitochondria, chloroplasts, and the cytoplasm.

3.2. Conserved Domain and Phylogenetic Analysis of Hevea bHLH

Analysis of conserved motifs within the HbbHLH family revealed a notable level of conservation, with 23 amino acid residues exhibiting over 50% conservation and 10 residues showing more than 75% conservation across all proteins (Figure 2). The majority of the conserved residues were consistent with previous research, particularly in regions associated with DNA-binding and dimerization. For instance, in the basic region of animal bHLH proteins, all bHLH proteins binding to the E-box (5'-CANNTG-3') have glutamic acid at position 9[53]. In rubber trees, the critical E9 residue is highly conserved and is found in 86% of the analysed proteins; 62% of rubber tree bHLH proteins possess the characteristic configuration H5-E9-R13 of Group B animal bHLH proteins[6], which has been shown to bind to the variant of the E-box (5'-CACGTG-3')[54]. Highly conserved hydrophobic residues in helices 1 and 2 are essential for dimerization[55]. In rubber trees, leucine (L) residues are present at position 23 in each bHLH protein, suggesting that this residue may play a significant role in dimerization.

To classify the bHLH gene family of the rubber tree and identify its evolutionary relationship with the Arabidopsis bHLH family, a phylogenetic tree was constructed using 180 rubber tree bHLH protein sequences and 121 Arabidopsis bHLH protein sequences (Figure 3). According to the classification criteria of the AtbHLH family by Heim et al.[6], 180 rubber tree bHLH family members were divided into 23 subfamilies. Subfamily XII was the largest, containing 26 HbbHLH proteins, while the smallest subfamilies, VIIIa and IIIe, contained only two HbbHLH proteins. No rubber tree bHLH proteins were classified into the families IVd and VI, suggesting that these genes were likely lost during evolution.

3.3. Gene Structural, Regulatory Elements, Chromosomal Location, and Synteny Analysis of Hevea bHLH

The gene structure of the Hevea bHLH family exhibited substantial variation, with intron numbers ranging from 0 to 12 across the 180 identified genes (Figure S1). Notably, 20 genes are intronless, primarily distributed across subfamilies IIId, IIIe, and VIIIb, suggesting potential functional specialization or evolutionary pressure favoring rapid transcriptional responses[56]. Members of the same subfamily have similar structures. For example, all members of the IIId family, except HbbHLH109, were intronless (Figure 4a). Families Ia, Ib, IIIa, IIIb, IVc, IVb, Vb, VIIIa, and X contain 1-4 introns, while subfamilies such as IIIf, Va, VIIa, VIIb, and XII exhibit more complex structures with 5 to 10 introns (Figure S1). Conserved motif analysis using the MEME tool revealed that all HbbHLH genes share a core motif (motif 3), most of genes contain 3-5 conserved motifs, while members of families IIId, IIIe, and IIIf contain more conserved motifs, ranging from 7 to 11. The most conserved motifs are distributed at the C-terminus, such as those in families I, IIIb, IVa, VIIb, VIIa, XII, XI, IX, VIIIa, VIIIb, and VIIIc. Families IIId, IIIe, and IIIf had a segment of conserved motifs at the N-terminus (Figure S1).

An analysis of cis-regulatory elements within the 1500 bp promoter regions of the HbbHLH genes uncovered a rich diversity of regulatory motifs (Figure S2), reflecting the complex regulatory networks governing their expression. The upstream cis-regulatory elements of HbbHLH genes were divided into three main categories. The first category is related to abiotic stress, such as light-responsive elements (Box 4, G-box, GT1-motif), antioxidant response elements (ARE), low-temperature response elements (LTR), drought-induced elements (MBS), and wound response elements (WUN-motif, W box). The second category is related to phytohormone responsiveness, such as ethylene-responsive elements (ERE), abscisic acid-responsive elements (ABRE), methyl jasmonate-responsive elements (TGACG-motif), salicylic acid-responsive elements (TCA-element, SARE), gibberellin-responsive elements (P-box, GARE-motif), and auxin-responsive elements (TGA-element, AuxRR-core). The third category is related to plant growth and development, such as circadian rhythm regulatory elements (circadian), cell cycle regulatory elements (re2f-1), seed-specific regulatory elements (RY-element), endosperm-specific regulatory elements (AACA_motif), and meristem expression-related elements (CAT-box).

Among the 4187 cis-regulatory elements found, the first category had the highest frequency, with 2858 elements, accounting for 68% of the total. Light-responsive elements were particularly abundant, accounting for 46% of the total elements. The second category contained 1225 elements, mainly ethylene, abscisic acid, and methyl jasmonate-responsive elements. The third category of regulatory elements was less frequent with only 104 elements (Figure 4b). This suggests that HbbHLH genes may be primarily involved in the responses to various plant hormones and abiotic stresses.

Chromosomal localization results showed that the 180 HbbHLH genes were unevenly distributed across 18 chromosomes of the rubber tree (Figure 4c). Most of the HbbHLH genes were concentrated in the terminal regions of the chromosomes; chromosomes 10 and 15 had the fewest, with only four genes, while chromosome 5 had the most, containing 23 HbbHLH genes. Among the 180 HbbHLH genes, 12 tandemly duplicated genes formed five gene clusters, which are marked in green in Figure 4c.

Synteny analysis indicated that segmental duplications were the predominant factor driving the expansion of the bHLH gene family in the rubber tree. Among the 180 identified HbbHLH genes, 136 (approximately 75%) were segmentally duplicated, underscoring the critical role of genome-wide duplication events in shaping the diversity of this transcription factor family (Figure 5). These results suggest that segmental duplications significantly contributed to the retention and diversification of HbbHLH genes, allowing for functional divergence and adaptation in response to the unique selective pressures experienced by rubber trees.

3.4. Gene Ontology and Interaction Network Analysis of Hevea bHLH

Based on previous analysis of HbbHLH gene structure and motif diversity, as well as literature mining, it is speculated that HbbHLH genes are involved in a range of biological processes. To further investigate the function of these genes, Gene Ontology (GO) enrichment analysis was performed. The results showed that 94 HbbHLH genes were enriched among 134 GO terms. Specifically, all nine molecular function (MF) terms were associated with DNA binding and RNA polymerase II specificity, whereas the two cellular component (CC) terms were related to transcription regulation complexes. More importantly, enrichment in 123 biological processes (BP) highlights the broad functional repertoire of HbbHLH genes, with notable processes including positive regulation of RNA polymerase II-mediated transcription, plant development and differentiation, light response, and morphogenesis (Figure 6). These findings suggest that HbbHLH genes are integral to the coordination of growth and development of rubber trees.

Furthermore, various biological processes, including reproductive development, seed and seedling development, hormone biosynthesis (for example brassinosteroids and gibberellins), metabolism, secondary metabolism (mainly flavonoids), pigment metabolism, iron ion homeostasis, metabolism, and photosynthesis regulation were enriched (Table S2). These findings indicate that bHLH transcription factors are broadly involved in critical biological processes in rubber tree reproduction, growth, development, and stress responses, highlighting their significant roles.

In A. thaliana, the functions and interactions of various bHLH proteins have been well characterized by mapping Hevea bHLH proteins onto their Arabidopsis homologues, and an interaction network for 65 HbbHLH genes has been inferred. STRING analysis revealed extensive interactions among these HbbHLH proteins (Table S3), consistent with the characteristic ability of the bHLH family to form homodimers or heterodimers[57]. For example, homologues of SCRM (HbbHLH004), MUTE (HbbHLH126), FAMA (HbbHLH062), and SPCH (HbbHLH049) interact to regulate stomatal development[58]. Additionally, IDT (HbbHLH063), PYE (HbbHLH038), and ILR3 (HbbHLH061) interact with each other and contribute to iron homeostasis[59]. These findings provide valuable insights into the potential functions of the bHLH genes in rubber trees.

The predicted interactions between HbbHLH proteins and other genes revealed three major functional clusters: one involving iron metabolism, with interactions including BTS, FRO2, ZIF1, and NAS4; another related to plant growth and development, particularly brassinosteroid signalling (for example, IBH1, BZR1), and light responses (for example, PAR1 and CRY2); and a third cluster associated with epidermal and stomatal development, involving genes such as EPFL4 and EPF2 (Figure 7). These interaction patterns suggest that the bHLH gene family plays a significant regulatory role in metal homeostasis, hormone signalling pathways, light-responsive processes, and epidermal and stomatal development in rubber trees.

3.5. Expression Patterns of Hevea bHLH in Diverse Environmental and Physiological Events

To investigate the spatiotemporal expression patterns of HbbHLH genes, expression levels (FPKM) from different tissues and treatments of 16 rubber tree varieties were downloaded from the HeveaDB database[51] (Table S4). 39 HbbHLH genes with expression levels below 2 in all samples, such as HbbHLH009, HbbHLH010, and HbbHLH019, were filtered out and a heatmap was created (Figure 8a).

Most HbbHLH genes exhibit significant tissue-specific expression in various tissues. For example, HbbHLH015 and HbbHLH075 showed significantly higher expression in latex than in other tissues; HbbHLH060, HbbHLH170, and HbbHLH108 were highly expressed in male flowers, and HbbHLH167, HbbHLH123, and HbbHLH096 were highly expressed in the bark. Genes within the same subfamily displayed similar expression patterns. Except for HbbHLH177 and HbbHLH178, genes in the IIId and IIIe families were significantly expressed in latex. Most genes highly expressed in male flowers belonged to the XII family, whereas those highly expressed in the bark were mainly from the Ib family.

In rubber trees, jasmonic acid (JA) and ethylene are crucial for promoting laticifer differentiation[60] and latex production[61]. Previous studies have indicated that certain bHLH transcription factors respond to JA signals and integrate into an ethylene-regulated signalling network. In this study, we examined the expression patterns of HbbHLH genes under JA and ethylene stimulation using transcriptome data. The analysis revealed that most HbbHLH genes that were highly expressed in latex were sensitive to JA and ethylene. For example, HbbHLH015, HbbHLH074, and HbbHLH164 were significantly upregulated by ethylene and JA treatment. Additionally, as treatment duration increased, these HbbHLH genes exhibited specific expression trends, revealing the dynamic characteristics of gene expression under hormone regulation.

To further understand the dynamic changes in gene expression during ethylene treatment, Mfuzz was used to group all differentially expressed genes at different time points into eight clusters (Figure 8b). Clusters 1 and 2 showed continuous upregulation and downregulation, respectively, after 12 h of ethylene treatment. Changes in Clusters 3 and 4 occurred at approximately 3 h, with significant downregulation and upregulation observed afterwards, respectively. Clusters 5 and 6 exhibited significant upregulation, followed by downregulation at approximately 12 h and 3 h, respectively. Cluster 7 displayed a general downregulation trend before 12 h and significant upregulation afterwards. Cluster 8 was significantly upregulated before 3 hours of treatment, followed by fluctuations between downregulation and upregulation.

Given the significance of cold resistance mechanisms in rubber tree research and the role of bHLH transcription factors in response to cold, a similar time series analysis was conducted on HbbHLH genes under cold treatment. Mfuzz divided them into six clusters (Figure 8c), revealing that the genes in clusters 1 and 4 were significantly downregulated and upregulated after 8 h, respectively. Genes in cluster 2 were consistently downregulated throughout the cold stress period; cluster 3 showed significant downregulation between 0 and 2 h; cluster 6 showed significant upregulation between 0 and 2 h, followed by significant downregulation; and cluster 5 exhibited an overall upregulation trend.

Additionally, HbbHLH genes exhibited significant expression dynamics at various stages of leaf growth and development in rubber trees. Specifically, HbbHLH077 exhibited the highest expression during Stage B, whereas HbbHLH171 peaked during Stage BC. The expression levels of HbbHLH083, HbbHLH076, HbbHLH006, etc. gradually increased with leaf maturation. Conversely, HbbHLH045 and HbbHLH133 showed decreasing trends throughout the developmental stages. Notably, most HbbHLH genes were abundantly expressed during Stages B, BC, and C but significantly decreased at Stage D; Furthermore, members of the HbbHLH family also responded to the disease. During tapping panel dryness (TPD), a common stress condition encountered by rubber trees due to frequent tapping for latex collection, genes such as HbbHLH156, HbbHLH176, and HbbHLH086 were upregulated, whereas HbbHLH096 and HbbHLH167 were downregulated. When confronted with Brown Blast infection, the expression of genes such as HbbHLH094 and HbbHLH163 decreased, but HbbHLH059, HbbHLH123, HbbHLH133, and others were upregulated, while under drought stress, the expression of genes such as HbbHLH057, HbbHLH164, and HbbHLH167 was inhibited, while genes such as HbbHLH029, HbbHLH051, and HbbHLH176 were upregulated. These observations strongly support the functional diversity of bHLH transcription factors, indicating their extensive involvement in regulating plant growth and development, hormone responses, and stress adaptation.

Genes that were significantly highly expressed in all tissues or showed differential expression in at least one of the five tissues were selected for RT-qPCR validation of the expression levels in different tissues (bark, leaves, latex, female flowers, and male flowers). The qPCR results were largely consistent with the transcriptome data. HbbHLH130 showed a relatively uniform expression across all tissues. HbbHLH015 was highly expressed in the latex. HbbHLH060 was highly expressed in the male flowers. HbbHLH083 was significantly highly expressed in the leaves and almost undetectable in other tissues. HbbHLH123 showed higher expression levels in bark and flowers (Figure 9).

Based on the transcriptome data and previous research, 11 genes were selected for RT-qPCR to further elucidate the expression patterns of HbbHLH genes under ethylene and methyl jasmonate treatments. Some of these selected gene homologues have been shown in other plants to be involved in responses to ethylene and methyl jasmonate.

The results indicated that most of the selected genes responded to ethylene and methyl jasmonate treatment. Under ethylene treatment, HbbHLH015, HbbHLH115, HbbHLH163, and HbbHLH164 displayed an initial upregulation, followed by downregulation, with a peak at 8 h post-treatment. HbbHLH098, HbbHLH109, and HbbHLH116 were significantly upregulated within the first 4 h, followed by downregulation, with expression levels rebounding 12 h post-treatment. HbbHLH059, HbbHLH074, HbbHLH075, and HbbHLH094 exhibited fluctuations in expression during the early stages of treatment, peaking at 24 h (Figure 10a). Under methyl jasmonate treatment, most genes were significantly upregulated, peaking at 8-12 hours. Notably, HbbHLH094, HbbHLH163, HbbHLH074 and HbbHLH115, demonstrated similar expression patterns under JA treatment (Figure 10b). Additionally, HbbHLH075 and HbbHLH163 showed similar transcriptional changes after treatment with ethylene and methyl jasmonate. However, the quantitative results for individual genes did not always align with the transcriptional changes observed in transcriptome data. This discrepancy may be due to sample heterogeneity, experimental variability, or differences in the data normalization methods.

Similarly, based on transcriptome data and previous research, 12 HbbHLH genes were selected to validate changes in their expression patterns under cold conditions. The results showed that genes such as HbbHLH094, HbbHLH076, and HbbHLH163 were significantly upregulated under cold conditions, whereas HbbHLH043, HbbHLH051, and HbbHLH179 were significantly downregulated. Interestingly, HbbHLH094 and HbbHLH163, which exhibited similar expression trends under methyl jasmonate treatment, also displayed similar expression patterns under cold conditions (Figure 11).

4. Discussion

4.1. Evolutionary Analysis of the bHLH Family

In most species, despite varying selective pressures, the number of bHLH genes remains relatively consistent. For example, the 119.1 Mb genome of A. thaliana has 162 bHLH genes[7]. In contrast, the 385.7 Mb genome of rice (Oryza sativa) contains 167 genes[9]. With a 392.2 Mb genome, Populus trichocarpa has 183 bHLH genes[8]. A total of 197 of these genes were found in Pyrus × bretschneideri[62], which had a genome size of 509.1 Mb. Finally, H. brasiliensis has 180 bHLH genes despite having a larger genome (1.9 Gb). This stability suggests that bHLH transcription factors play essential roles in fundamental biological processes, and that the family has undergone considerable evolutionary conservation.

In the rubber tree, the bHLH genes were categorized into 23 subfamilies, supported by phylogenetic analysis, gene structure, and conserved motifs. The number of members in most subfamilies was proportional to that in Arabidopsis, and most subfamilies contained Arabidopsis members. This cross-species consistency implies that bHLH genes have evolved conservatively and their functions may have been preserved over evolutionary time. Consequently, members of the same subfamily, from different species, are likely to perform similar functions.

4.2. Functional Prediction of Hevea bHLH Genes

The bHLH transcription factors have been characterized in various plant species, and their diverse functions include plant growth and development, hormone responses, and stress resistance. In H. brasiliensis, we identified bHLH genes with conserved structural domains and diverse biological process clusters, suggesting that these genes also exhibit a broad functionality in rubber trees. Genes with close phylogenetic relationships and protein domains often have similar functions[63]. Previous studies have shown that in Arabidopsis, three bHLH proteins from the Ia family, SPCH, MUTE, and FAMA, along with their heterodimeric partners SCRM, collectively mediate cell-state transitions and induce stomatal development[64,65]. In rubber trees, homologous genes HbbHLH049, HbbHLH126, HbbHLH062, and HbbHLH004 are likely to have similar functions. This conclusion is supported by GO enrichment analysis, which showed that these four genes, along with other members of the Ia family, are enriched in the biological process of stomatal complex development. Similarly, the IX subfamily in Arabidopsis includes AtbHLH80 (FBH) genes, which regulate flowering time via the CO pathway[15]. In Hevea, the IX subfamily members HbbHLH076, HbbHLH118, and HbbHLH139 were enriched in flowering-related GO terms (GO:0048573), suggesting similar regulatory roles. Additionally, the Ib subfamily members HbbHLH038, HbbHLH120, and HbbHLH137 may regulate iron homeostasis in a manner similar to AtbHLH38 (ORG2) and AtbHLH39 (OGR3)[18].

Protein-protein interaction analysis identified interactions between bHLH genes and other functionally characterized genes. For example, HbbHLH156 interacts with 16 other bHLH proteins across various subfamilies, suggesting that it is involved in important regulatory networks. However, research on this gene remains limited and further studies are warranted to explore its role in rubber tree physiology.

4.3. Potential Role of IIId and IIIe Subfamilies in Latex Biosynthesis in Hevea

Phylogenetic analysis revealed an expansion of the IIId subfamily in H. brasiliensis, with four times more members than those in A. thaliana. Nearly all members of this subfamily are highly expressed in latex, particularly HbbHLH015 and HbbHLH075, which show over tenfold higher expression in latex than in other tissues. Within the IIIe subfamily, HbbHLH094 and HbbHLH059, homologues of Arabidopsis MYC genes, were also highly expressed, specifically in latex. In Arabidopsis, MYC genes regulate jasmonic acid (JA)-induced processes, and in Hevea, MYC homologues are involved in rubber biosynthesis via the COI1–JAZ–MYC2 signaling pathway[66]. Further analysis of the transcriptional responses of these subfamilies to JA and ethylene treatments revealed that most members responded to both the treatments. In particular, significant transcriptional changes were observed under jasmonic acid treatment, where HbbHLH163 and HbbHLH094, as well as HbbHLH059, HbbHLH074, and HbbHLH115, exhibited similar expression patterns. This suggests that these genes participate in the JA-mediated regulation of rubber biosynthesis.

Enrichment analysis also linked the genes in group III to the synthesis and metabolism of secondary metabolites, such as anthocyanins and flavonoids in subfamilies IIIf and IIIe and lipids in subfamily IIIb. Given the close association between natural rubber production and secondary metabolites, particularly terpenoids, it is possible that IIId subfamily members are involved in latex production through their roles in secondary metabolite synthesis and metabolism. Further investigations are required to elucidate the specific functions of these subfamilies in rubber biosynthesis.

4.4. Regulatory Mechanism of bHLH Response to Cold Stress

Several HbbHLH genes, including HbbHLH043, HbbHLH051, HbbHLH076, HbbHLH094, HbbHLH116, and HbbHLH163, showed significant transcriptional changes after cold treatment, although no cold stress-related GO terms were enriched. According to the GO enrichment analysis, these genes are involved in multiple subfamilies and are implicated in a variety of biological processes. This suggests that these genes may participate in indirect pathways related to cold stress, such as light signalling, stomatal movement, or hormone response pathways.

Previous studies have shown that MYC2 mediates cold stress responses through the JA signalling pathway in apples[67]. In Hevea, HbbHLH094 exhibited significant transcriptional changes under both low-temperature and JA treatments, suggesting that it might also participate in JA-mediated cold stress responses. In contrast, the segment duplication gene HbbHLH059 did not display similar changes in expression, likely because of differences in cis-regulatory elements. Notably, HbbHLH094 and HbbHLH163 displayed similar transcriptional patterns in response to both methyl jasmonate and cold treatments, indicating that they may jointly participate in JA pathway regulation of cold stress responses.

5. Conclusions

This study identified 180 bHLH genes in H. brasiliensis, thereby highlighting the evolutionary conservation and functional diversity of the bHLH gene family. Through the integration of extensive transcriptomic data, GO enrichment analysis, and protein interaction analysis, we demonstrated that bHLH genes in Hevea are implicated in a broad spectrum of biological processes, including growth and development, photomorphogenesis, iron ion homeostasis, secondary metabolite synthesis and metabolism, and responses to abiotic stress and hormonal signals. Our analysis identified two subfamilies, IIId and IIIe, that are potentially crucial for natural rubber biosynthesis. Furthermore, several bHLH genes exhibited significant transcriptional changes in response to cold stress, suggesting their involvement in stress-adaptation mechanisms. This study provides guidance for the functional characterization of individual family members and offers potential targets for genetic manipulation aimed at enhancing stress tolerance and improving the growth performance of rubber trees.