Genomic Survey and Cold-Induced Expression Patterns of bHLH Transcription Factors in Liriodendron chinense (Hemsl) Sarg.

Abstract

bHLH transcription factors play an animated role in the plant kingdom during growth and development, and responses to various abiotic stress. In this current study, we conducted, the genome-wide survey of bHLH transcription factors in Liriodendron chinense (Hemsl) Sarg., 91 LcbHLH family members were identified. Identified LcbHLH gene family members were grouped into 19 different subfamilies based on the conserved motifs and phylogenetic analysis. Our results showed that LcbHLH genes clustered in the same subfamily exhibited a similar conservative exon-intron pattern. Hydrophilicity value analysis showed that all LcbHLH proteins were hydrophilic. The Molecular weight (Mw) of LcbHLH proteins ranged from 10.19 kD (LcbHLH15) to 88.40 kD (LcbHLH50). A greater proportion, ~63%, of LcbHLH proteins had a theoretical isoelectric point (pI) less than seven. Additional analysis on the collinear relationships within species and among dissimilar species illustrated that tandem and fragment duplication are the foremost factors of amplification of this family in the evolution process, and they are all purified and selected. RNA-seq and real-time quantitative PCR analysis of LcbHLH members showed that the expression of LcbHLH35, 55, and 86 are up-regulated, and the expression of LcbHLH9, 20, 39, 54, 56, and 69 is down-regulated during cold stress treatments while the expression of LcbHLH24 was up-regulated in the short term and then later down-regulated. From our results, we concluded that LcbHLH genes might participate in cold-responsive processes of L. chinense. These findings provide the basic information of bHLH gene in L. chinense and their regulatory roles in plant development and cold stress response.

1. Introduction

Globally out of all abiotic stress factors, cold, drought, and heat stresses are declared as the most complex ones affecting plant growth, survival, and crop productivity. Molecular regulation at the post-transcriptional level possesses a vital role for development, growth, nutrient allocation, and defensive mechanism in plants [1,2]. The bHLH family regulates growth and development, morphogenesis, and stress responses in plants [3,4,5], characterized by a helix-loop-helix (HLH) domain, with an approximated 15 amino acids N-terminal as the base region: known for recognizing and binding to specific DNA while, the C-terminal is the HLH region with about 50 amino acids [6,7,8]. The helix is also associated with DNA sequences that recognize protein-specific binding [9] and can form homodimer or heterodimer with other proteins [10]. On top of an α-helix near the N-terminal is another α- helix [11]. The two α-helices are connected with a ring formed by amino acid chains to form an HLH structure.

Generally, bHLH transcription factors are known to act as transcriptional activators or inhibitors for seed germination and flowering regulation [9]. However, a study in Arabidopsis mutant srl2, AtPIF4 (AtbHLH09) spectacled a specific role in the signaling network in phytochrome B (phyB) and in light regulation [12]: AtPRE1 (AtbHLH136) and ILI1 were also identified to regulate cell elongation by interacting with IBH1 (AtbHLH158) under the action of brassinosteroids (BR) and gibberellin signals [13]. Moreover, AtPRE1 (AtbHLH136) and IBH1 (AtbHLH158) form regulatory system with AtACE1/2/3 (AtbHLH049/074/077) that competitively regulate cell growth. IBH1 (AtbHLH158) has also been shown to negatively regulate cell growth by interacting with the positive regulatory gene AtACE1/2/3 (AtbHLH049/074/077) [13]. Certain members of the bHLH transcription factor family have also been shown to enhance resistance to harsh conditions when plants retort to abiotic stresses [14,15]. For instance, overexpression of AtICE1 (AtbHLH116) and AtICE2 (AtbHLH33) can augment the expression of CBF promoter at low temperature and mend the stress resistance of transgenic plants [16,17]. Feng et al. [18] has also demonstrated that MdCIbHLH1 protein binds to the MdCBF2 promoter and upregulates the expression of CBF2 through the C-repeat-binding factor (CBF) pathway and promote the cold tolerance of transgenic apple plants. A study in trifoliate orange has also shown PtrbHLH to increase cold resistance by activating PtrCAT [19].

To date, research on different plant genomes has concurred that the bHLH transcription factor family is incessantly distinguished, with the structural characteristics and response profiles to various environmental stresses [10,20,21,22]. Nonetheless, few studies on the bHLH gene family of the forest tree species have been conducted with less on the L. chinense. L. chinense is a kind of tall deciduous tree, which is of economic, ornamental, medicinal, and ecological value [23,24]. The recent release of the L. chinense genome provided the opportunity for its LcbHLH gene family (which will be referred to as Lc in this study) to be analyzed [23]. In this current study we identified 91 LcbHLH transcription factors, which were further analysed using Bioinformatic approach for evolution, conserved motif arrangement, exon-intron patterns, and other physiochemical proprieties. Additionally, each subfamily of the LcbHLH gene family was shown to play imperative biological functions in abiotic stress responses. The identification and distinctive analysis of the bHLH transcription factor of L. chinense will assist in comprehending the structural characteristics of gene families in L. chinense and preliminarily predict the function of bHLH members, which will provide the gene resources for the improvement of L. chinense germplasm by genetic engineering technology in the future.

2. Materials and Methods

2.1. Identification and Physicochemical Properties Analysis of bHLH Family Members of Liriodendron chinense

The nucleic acid and protein sequences of L. chinense were collected from the local protein database [23]. The protein sequences of the bHLH family of Arabidopsis and rice were retrieved and downloaded from the plant transcription factor database (http://planttfdb.cbi.pku.edu.cn (accessed on 12 November 2021)) [25]. The bHLH protein sequences of Arabidopsis and rice were used as query sequences, while the candidate protein-containing bHLH/HLH domain was screened from the L. chinense database by local blastp program. Then, the HMMER model downloaded from the Pfam database was used to identify the candidate bHLH protein of L. chinense in a local protein database. Finally, proteins with the bHLH/HLH domain were taken as the final bHLH family members of L. chinense. The physical and chemical properties (including molecular weight, isoelectric point, and hydrophilicity) of LcbHLH family members were analyzed using the Protparamin EXPASY database.

2.2. Phylogenetic Analysis of LcbHLHs

ClustalX2 was used for multiple sequence alignment of the bHLH domain. The bHLH proteins of three plants, rice, Arabidopsis, and poplar, have been downloaded from National Center for Biotechnology Information (NCBI). The phylogenetic tree was constructed using MEGA7.0 with the Neighbor-Joining method [26,27]. The evolutionary distance was obtained through the p-distance method, with the distances employed to estimate the number of amino acids at each locus. The reliability of each phylogenetic tree was guaranteed by 1000 bootstrap sampling iterations.

2.3. Chromosome Location and Gene Replication of LcbHLHs

The data of the chromosomal location of LcbHLH members were obtained from annotated files in the Liriodendron genomic database, while the distribution of LcbHLH members was plotted using the biological software TBtools [28]. The gene replication events were analyzed according to the following three standard definitions: (1) the length of one shorter sequence is greater than 70% of that of the other longer sequence; (2) the similarity between the two sequences is greater than 70%; (3) two genes separated by five or fewer genes in a 100 kb chromosome segment are considered as tandem repeat genes [29]. To analyze the collinearity correlation between LcbHLHs and bHLHs in other species, the genome data of Arabidopsis and rice were downloaded from Ensemble (http://plants.ensembl.org/index.html (accessed on 13 November 2021)). The multicollinearity scanning tool MCsanX was employed to compare the whole genome sequence of Liriodendron with that of Arabidopsis and rice, respectively [30]. The visualization of chromosome distribution was obtained through the Circos in TBtools. The ratio of Ka/Ks was calculated by using KaKs_calculator to acquire the natural purification selection between target gene pairs [31].

2.4. Analysis of Gene Structure, Conserved Motifs and Cis-Regulation Elements of LcbHLHs

TBtools software was adopted to map the gene structure of LcbHLH members onto a diagram. MEME was used to predict and analyze the conservative motif of the bHLH protein in L. chinense. Cis-regulation elements of LcbHLH members were predicted by the software Plantcare (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/ (accessed on 25 November 2021)) and plotted by TBtools.

2.5. Analysis of Protein Interaction among LcbHLHs

The protein interaction network was generated using the STRING (www.string-db.org (accessed on 3 December 2021)) based on the high homology between LcbHLHs and AtbHLHs proteins. In addition, six LcbHLH proteins with high homology to AtbHLH were selected to map the extrafamilial protein interaction network using Cytoscape 3.8.2 [32].

2.6. Three-Dimensional Structure Modeling and Verification of bHLH Protein

The full-length atomic structures of LcbHLH24, LcbHLH72, and AtICE1 proteins were constructed based on the synthesis method on the Robetta online website. Homologous modeling was used for proteins with the sequence matching model, while the threading method was used for proteins with the sequence non-matching model. Then, the sequence was assembled to construct the protein structure. The reliability of their protein structures was further confirmed by ERRAT, PROVE, and Ramachandran on the online website Savesv6.0. VMD software was used for 3D modeling.

2.7. Expression Analysis of LcbHLHs in Response to Cold Stress by RNA-seq and qRT-PCR

The somatic embryo-regenerated seedlings of hybrid Liriodendron with consistent growth were cultured in an incubator (23 °C, 16 h light, and 8 h dark) and then treated at 4 °C. Seedling leaves were sampled at 0 h, 6 h, 1 day, and 3 days with three biological replicates. The collected leaves were quickly frozen in the liquid nitrogen and put in a −80 °C refrigerator for storage. Transcriptome sequencing was performed on the above samples. Transcriptome data of LcbHLH members were extracted from the sequencing results. The expression levels of each member at each period of cold stress treatment (the maximum expression value of each LcbHLH gene was set to 1, and then the expression values of the gene at other stress and growth stages were normalized to the maximum expression value) were normalized and displayed on the heatmap. The expression patterns of ten LcbHLH members were determined by quantitative RT-PCR analysis (qRT-PCR). The qRT-PCR was performed using SYBR-green in the Roche Light Cycler^®480 real-time PCR system (Switzerland, Sweden). The relative expression abundance of LcbHLH was calculated with the ∆∆CT method. 18s rRNA was used as the internal reference. All qRT-PCR primers were designed by Primer5.0 and were listed in Table S1.

3. Results

3.1. Identification and Physiochemical Characteristics of LcbHLHs

Based on the search of the conserved bHLH domain (Pfam number: PF00010), 91 LcbHLH family members were recognized after further validation in the conserved domain database (CDD) and Pfam database. They were renamed as LcbHLH 1~91 based on their chromosomal position. The physical and chemical properties of LcbHLH members were computed. Analysis of the hydrophilicity value of all LcbHLH proteins showed a negative total average value that ranged from −0.816 (LcbHLH47) to −0.143 (LcbHLH68), concluding that LcbHLH proteins are hydrophilic. The Molecular weight (Mw) of LcbHLH proteins ranged from 10.19 kD (LcbHLH15) to 88.40 kD (LcbHLH50), the majority (61%) were in the range of 21.41 kD to 48.85 kD, and the molecular weight of 24 members (about 26%) was in the range of 20 kD to 30 kD. Additionally, the theoretical isoelectric points (pI) of LcbHLH proteins ranged from 4.59 (LcbHLH81) to 9.91 (LcbHLH53). Most LcbHLH proteins (about 63%) were less than 7, and about 30% of LcbHLH proteins had a pI between 6 and 7 (Table S2).

3.2. Phylogenetic Characteristics of LcbHLHs

To fully comprehend the evolutionary relationship of the identified LcbHLH protein sequences, L. Chinense (Lc), Arabidopsis thaliana (At), Oryza sativa (Os), and Populus. trichocarpa (Pt), bHLH gene families were further compared and subjected in phylogenetic tree analysis (Figure 1A and Figure S1). A total of 581 bHLH protein sequences were obtained and divided into 31 groups, which were identified as evolutionary branches with high bootstrap values. Among the 31 subfamilies, 26 subfamilies were presented in all four species, signifying that the genes of these subfamilies had high homology in the four species and strong phylogenetic conservatism. Some LcbHLH genes in Arabidopsis and rice were clustered in the same subfamily. LcbHLH proteins were clustered in 29 subfamilies and an orphan sequence was observed. Subfamily 13 was clustered with LcbHLH14, LcbHLH15, AtPRE1/2/3/4/5, and AtKDR. Subfamily 17 was clustered with LcbHLH24, LcbHLH31, LcbHLH82, LcbHLH18, At033SCRM, and At116ICE1 (Table S3). Subfamily 25 was clustered with LcbHLH16, LcbHLH69, LcbHLH78, AtSPCH, and OsSPC1/2. Additionally, subfamily 6 was only found in Poplar, indicating individual evolution and functional diversity of Poplar (Figure 1B).

3.3. Gene Structure and Conserved Motifs of LcbHLHs

Gene structure prediction plays an animated role in studying the evolution of gene family members. To further explore the phylogenetic relationships within the LcbHLH members, the intron/exon structures of the LcbHLH gene were analyzed based on the genomic annotation files of 91 LcbHLH members in combination with phylogenetic tree (Figure 2A). The number of introns in the LcbHLH gene ranged from 1 to 11. LcbHLH genes were clustered together by parallel exon/intron patterns in exon length and intron number (Figure 2B).

In this study, the configuration of the LcbHLH conservative motif was discovered through the protein conservative theme sites predicted by online software MEME (Figure 1C and Table S4). bHLH conserved domain was constituted by motif 1 and motif 2 (Figure 1C). The meticulously connected LcbHLH proteins on immediate evolutionary branches of the phylogenetic tree had the same or comparable motif structures. Moreover, there were significant differences between dissimilar subfamilies, suggesting that members of the identical subfamily of bHLHs might play related roles in L. chinense. Seven subfamilies shared motif 11, eleven subfamilies shared motif 3, and nine subfamilies shared motif 4. Motif 19 only occurred in subfamily 4, motif 17 and motif 20 only occurred in subfamily 10, motif 16 only occurred in subfamily 11, motif 10 only occurred in subfamily 12.

3.4. Cis-Regulation Elements of LcbHLHs

The cis-regulatory element plays an imperative role in regulating the expression of stress response genes [33]. The presence of the cis-elements of the LcbHLH members in the promoter region (2000 bp upstream of the transcription initiation site) was predicted. Twenty-five typical elements with relatively robust functions were divided into three groups shown in Figure 3. Based on the functional annotations, cis-elements were categorized into three major classes: plant growth and development, phytohormone responsive, and abiotic and biotic stresses (Figure 3). Our findings showed that G-Box and ABRE were the most represented transcription factors in the LcbHLH gene family. Specifically, LcbHLH7 had the most representation of G-Box and ABRE. 67 LcbHLH members had elements responsive to the methyl Jasmonate, including CGTCA-motif and TGACG-motif. Fifty-four members had gibberellin-responsive elements, including P-box and GARE-motif. Seventy-two members had salicylic acid responsiveness elements, TCA-element. Moreover, 45 members had auxin-responsive elements, including AuxRR-core and TGA-element. 52 LcbHLHs contained LTR elements that might be interrelated to the cold stress response of L. chinense.

3.5. Intergenomic Collinearity and Gene Replication of LcbHLHs

Amongst 91 LcbHLH genes, 89 were distributed on 19 chromosomes, and the other two were assigned to unassembled genomic contigs (Figure 4). The number of LcbHLH genes on each chromosome ranged from 1 to 9.

The analysis of genome-wide replication, fragment replication, and tandem replication of gene family has a significant role in explaining the process of gene family expansion. In this analysis, intraspecies comparisons of L. chinense and A. thaliana, L. chinense, and rice were implemented at the genome-wide level (Figure 5). A total of 24 pairs of replication genes were found in the LcbHLH family, and 21 pairs of gene clusters with high similarity were institute in LcbHLHs (Figure 5A). For example, the protein sequences of LcbHLH88 and LcbHLH89 shared 99.23% resemblance. The similarity between LcbHLH63 and LcbHLH62 was 99.65%, respectively.

Additionally, tandem repeat genes comprised the same number of exons due to closely related imitation associations. The tandem repeat genes LcbHLH14 and LcbHLH15 and LcbHLH62 and LcbHLH63 had a similar two exon and intron-exon structure pattern. Likewise, LcbHLH84 and LcbHLH85 had a similar intron structure pattern. Remarkably, as revealed in Figure 5A, there were four pairs of fragment-repetitive genes: LcbHLH3, LcbHLH4 and LcbHLH36, and LcbHLH37; LcbHLH12, LcbHLH13 and LcbHLH27, and LcbHLH28; LcbHLH69, LcbHLH70 and LcbHLH80, and LcbHLH82; LcbHLH47, LcbHLH48 and LcbHLH59, and LcbHLH60. Together, these results show that the LcbHLH gene family was amplified by fragment replication and tandem replication of the LcbHLH genes.

The tandem repeated bHLH gene has a related gene structure, motif composition, and expression. The tandem repeated and intra-and inter-chromosome repeated regions of LcbHLH members were examined in the present study. Our results showed that greater than 38% (15 tandem and 22 fragment-repeat genes) of the LcbHLHs might have evolved from some genomic replication event. The substitution rate (Ka/Ks) between nonsynonymous and synonymous was an operative quantity of selection pressure after gene replication [34]. Consequently, the Ka/Ks of the LcbHLH repeat gene was premeditated (Table S5). For all tandem repeat pairs, the Ka/Ks values were well below one, which indicated that there were purification options during amplification. Besides, for gene pairs with fragment repeats, all Ka/Ks were less than one, indicating that there was strong purification selection pressure during evolution.

With genome-wide comparison and analysis of L. chinense, A. thaliana, and rice, it was established that most LcbHLHs were positively homologous in rice and A. thaliana (54% and 60%), respectively (Figure 5B,C, Tables S6 and S7). The Ka/Ks ratios of L. chinense to rice and A. thaliana were 0.175 and 0.186, respectively. These results indicate that bHLH gene pairs underwent strong purification selection and that there was a close correlation between them before. In brief, gene replication events, including tandem and fragment repeats, appeared to be essential for the expansion of the bHLH gene family in Liriodendron, as well as for the functional preservation and differentiation.

3.6. Protein Interaction Network of bHLHs

Diverse bHLH proteins bind to specific DNA and regulate the downstream target’s transcription by forming homodimer or heterodimer mediated by their α-helix near the N-terminal [10]. Hence, protein interaction analysis is essential to fully review the function of LcbHLH proteins (Figure 6). It can be speculated that LcbHLHs might have played a role in forming protein complexes and attempted to construct an interaction network of LcbHLHs. In this current study, the interaction network within the LcbHLH gene family was constructed based on the orthogonal analysis of AtbHLHs (Figure 6A and Table S8). The protein interaction network indicated that most LcbHLH proteins could interact with more than one bHLH protein. More than a quarter of LcbHLH proteins can interact with four or more other bHLH proteins. Numerous imperious interactions were predicted, such as how CIB1 (LcbHLH7) can participate in the regulation of flowering time [35]. ICE1 (LcbHLH24, 31) interacts with FMA (LcbHLH32), SPCH (LcbHLH78, 79, 16) and MUTE (LcbHLH53) to regulate stomatal diversity [34]. LRL1 (LcbHLH75) and RDH6 (LcbHLH8) can interact with RSL2 (LcbHLH85 and 86) and contribute to the regulation of root hair development. These protein interaction networks further ascertained that the LcbHLH genes exerted their diverse biological functions through interaction and coordination with other members.

3.7. Structural Modeling of LcbHLH Protein

The bHLH transcription factor family plays a vital role in plant response to abiotic stress by forming dimer and its helical structure [36]. ICE, one of the bHLH families, activates CBF via transcription and persuades its expression, playing a central role in cold response and signal transcription [16,37,38,39,40,41]. The amino acid sequence of LcbHLH24 in L. chinense is extremely homologous to that of ICE1 in A. thaliana. For that reason, this research predicted that these two protein structures, LcbHLH24 and LcbHLH72 (homologous gene of AtRSL2), interacted with LcbHLH24 through RGE1 in the protein network (Figure 7A). The structure of LcbHLH24 consisted of 14 α-helices and 19 loops (Figure 7A), and the model of LcbHLH72 had ten α-helices and eight loops (Figure 7B). The three-dimensional structure of AtICE1 protein consisted of 14 α-helices and 12 loops (Figure 7C).

In the 3D model of LcbHLH24, the structural model could be roughly divided into three regions, exposed as a, b, and c. LcbHLH72 could be divided into two regions, designated as a and b. Three structural regions could be found in AtICE1, in which region b was similar to the structure of the other two proteins. Nevertheless, region a of LcbHLH24 and AtICE1 is a little richer than that of LcbHLH72. According to the homology modeling of SWISS-MODEL and the prediction of the conserved domain of NCBI (CDD), region b is the bHLH conserved domain of three proteins. The conserved structural region b of LcbHLH24 and AtICE1 was predicted by SWISS-MODEL to have the domain characteristics of the MYC2 subfamily. Alternatively, region b of LcbHLH72 showed high consistency with MITF/CLEAR box structure. Interestingly, special structural region Berninger c was only identified in LcbHLH24 and AtICE1, and region c in LcbHLH24 was almost identical to AtICE1. In summary, comparative analysis of LcbHLH24 and AtICE1 protein sequences, region c is a highly conservative Zipper domain.

3.8. Cold Stress-Induced Expression Pattern of LcbHLHs

The expression patterns of LcbHLHs under cold stress in transcriptome data were analyzed (Figure 8) to understand the responses of LcbHLHs to cold stress, and 78 LcbHLH genes were examined to express in the seedling leaves of L. chinense. During the cold stress treatment, the expression patterns of LcbHLH members were coarsely defined by constant up-regulations and down-regulations (Figure 8). The expression patterns under the cold treatment of 20 LcbHLH genes (22.2%) showed a constant up-regulation trend, 15 LcbHLH genes (16.7%) were incessantly down-regulated; 28 of the total LcbHLHs (31.1%) were up-regulated and then subsequently down-regulated with the extension of cold treatment time, and only four genes (4%) showed the down and then increased trends.

To further verify the expression pattern of LcbHLHs under cold stress, ten LcbHLHs (LcbHLH9, 20, 24, 35, 39, 54, 55, 56, 69, 86) were chosen to quantify the expression abundance in L. chinense by qRT-PCR. As shown in Figure 8B, the expression trends of these ten genes were almost consistent with their transcriptomic patterns. Three LcbHLH genes (LcbHLH35, 55, 86) showed an up-regulation trend in response to cold stress, six LcbHLH genes (LcbHLH9, 20, 39, 54, 56, 69) displayed a down-regulation trend, and the expression profile of one LcbHLH gene (LcbHLH24) was up-regulated at 1d and then down-regulated at 3d.

4. Discussion

Given the significant character and diverse functions in biological processes, the bHLH transcription factors have attracted more and more attention in recent years [21,42,43]. In this current study, members of the bHLH family identified from the genome of L. chinense had analogous structural characteristics to those of other species, especially the bHLH domain. That was highly conservative with 19 amino acid residues, of which five were base regions, five were distributed in the first helix, one in the loop, and eight in the second helix [44]. However, typical conserved sites were found in the domain of the L. chinense bHLH gene family, like the AtbHLH families. This indicated that LcbHLHs might have DNA-binding activity like that of A. thaliana [45].

We constructed a phylogenetic tree to better understand the evolutionary relationship of bHLH gene families between different species, L. chinense, A. thaliana, rice, and poplar. Interestingly, genes with the same functions were clustered into the same clade. For example, LcbHLH78, LcbHLH79 and AT5G53210 (AtSPCH), Os02g15760 (OsSPCH2) and Os02g33450 (OsSPCH1) were clustered into subfamily 25. We used this evolutionary clustering on the same branch to speculate the functional importance of identified LcbHLHs. Previous studies in A. thaliana have shown that AtSPCH can regulate the formation of stomata together with AtMUTE and AtFAMA [46]. In rice, SPCH and MUTE have also been shown to exhibit the same functional importance in stomatal formation [47]. Hence, it is reasonable to speculate that LcbHLH78 and LcbHLH79 are imperative genes regulating the stomatal switch in L. chinense. Furthermore, the LcbHLH24 and LcbHLH31 were also clustered into the same subfamily (subfamily 17) as AT1G12860 (AtICE1), AT3G26744 (AtICE2), Os11G32101 (OsICE1), and Os01G0310 (OsICE2). AtICE1 and AtICE2 are the main transcription factors found in A. thaliana responding to low-temperature stress [17]. OsICE1 can be phosphorylated by OsMAPK3, thus enhancing the activation of OsbHLH to its target gene OsTPP1 in response to low-temperature stress [48]. So, it is reasonable to speculate that LcbHLH24 and LcbHLH31 are most likely to participate in the signal transduction of L. chinense in response to low-temperature stress.

Similarly, exon-intron patterns and similar conservative motif arrangements are consistent with the subfamily classification. It is known that genes with few or no introns have low levels of expression in plants [49]. However, a gene structure with compact exons may facilitate rapid expression in response to both endogenous and exogenous stimuli [50]. We observed that the exon structures of LcbHLH5 and LcbHLH35 were relatively tight, and they both belonged to the subfamily 29. According to the transcriptomic data, the expression of these two genes under low-temperature stress was increasing in response to an increase in the duration of treatment exposure.

Genomic replication events occur throughout plant evolution, often leading to the expansion of gene families [51,52]. Tandem and fragment gene replication events are two major replication patterns common in the evolution of angiosperms [34,53] and play an essential role in gene family extension [51,54]. In the present study, several distinct gene clusters of LcbHLHs were distributed in the different chromosomes. Therefore, gene duplication might be an important reason for the large number of LcbHLHs. Gene replication is a common phenomenon in many organisms, which can regulate gene expression, improve genetic and environmental adaptability, and serve as a steppingstone in the evolution of new biological functions [55,56]. The relatively strong sequence diversity besides the bHLH domain suggests that the bHLH family has undergone extensive domain reorganization after gene replication [57]. More than 20 different conserved motifs with different arrangements were found in the bHLH family of L. chinense. Thus, extensive domain reorganization occurred in the protein structure of the bHLH members. This phenomenon implies that the evolutionary position of Liriodendron is difficult to determine accurately [23].

Time-specific expression patterns of genes in plant growth usually reflect variances in biological functions of gene family members and interactions among related pathways [58,59]. In transcriptional expression profiles, the diverse expression patterns of LcbHLH genes under cold stress inferred that each LcbHLH member might participate in the various cascades of signal transduction in L. chinense in response to cold stress. By predicting the cis-regulation elements of these LcbHLH genes, we observed regulatory elements responsive to temperature stress, including LTR, TCA, and AT-rich. The low-temperature responsive element LTR, with CCGAC as the core sequence, demonstrated diverse expression patterns under low-temperature stress, suggesting that LTR plays a key role in responding to low-temperature stresses [60,61]. CRT/DRE element is an important low-temperature response element in the bHLH family. CBF transcription factor can bind to CRT/DRE sequence and induce the expression of the COR gene to improve the cold resistance of plants [62,63].

Numerous proteins in the bHLH family are intricate in the tolerance to low-temperature stress, and ICE1 is a typical transcription factor that can regulate cold-responsive signal transduction in plants [37,64]. Two members (LcbHLH24 and LcbHLH31) were found to be highly homologous to AtICE1 and AtICE2 in the genome of L. chinense. The expression of LcbHLH31 was continuously up-regulated under cold stress, while the expression LcbHLH24 was continuously increased during one day but decreased after three-day treatment, but its abundance was still higher than that of the control. This indicated that two genes, LcbHLH31, and LcbHLH24, participated in the response of L. chinense to low-temperature stress. Over the comparative analysis of the protein sequences of LcbHLH24 and AtICE1, it can be inferred that LcbHLH24 has the characteristics of the typical ICE gene family, which contains an S-rich region and disulfide bonds. They can preserve the stability of its gene, but not in LcbHLH72. Consequently, it can be reasonably inferred that the stability of LcbHLH24 protein is stronger than that of LcbHLH72. Region c, which is found in the structure of LcbHLH24, shares the same characteristic with the structure of Zipper found in ICE of A. thaliana and other species. It can be expected and assumed that the special zipper protein structure of LcbHLH24 may be beneficial for further exploring and analysing the response of the bHLH family to low-temperature stress in L. chinense.

Protein-protein interaction analysis predicted interacted relationship among LcbHLHs, which of them were confirmed by previous reports. ICE1 [16], ICE2 [65], and MYB15 [66] have been recognized as regulatory factors that induce CBF expression. In response to low temperature, ICE1 can be sumoylated by SIZ1, thus promoting the binding of ICE1 and increasing CBFs expression [67]. In addition, SCRM2 plays an important role in regulating the stomatal development of SPCH, MUTE, and FAMA [36]. Evidence suggests that there may be a relationship between transcriptional regulation of environmental adaptation and stomatal development in plants [68].

5. Conclusions

This comprehensive genome-wide study systematically identified and functionally analyzed the bHLH gene family in L. chinense. A total of 91 LcbHLH family members were identified and divided into 31 subfamilies, which were unevenly distributed on 19 chromosomes of L. chinense. The reported gene structures, conservative motifs, and phylogeny further supported the characteristics of the phylogenetic trees. The amplification of the LcbHLH gene was due to duplication during evolution, suggesting that this gene family may play an important role in polyploid plants. Cis-regulation elements responding to low temperature were found in the upstream region of the LcbHLH gene, which indicated that the LcbHLHs might play an important role in response to cold stress. RNA-seq and qRT-PCR analysis showed that members of the LcbHLH genes had various expression patterns during cold treatments. These results may contribute to further functional studies of LcbHLH genes and may provide gene resources for the genetic improvement of L. chinense.