Genome-wide identification, molecular evolution, and expression analysis of the bHLH gene family in Mesembryanthemum crystallinum

Mesembryanthemum crystallinum (2n = 18), commonly known as ice plant, belongs to the family Aizoaceae. Native to southern and eastern Africa, M. crystallinum is now cultivated globally, and possesses high nutritional value for humans due to its rich content of antioxidants, such as phenolic compounds^[1]. Studies indicate its nutritional components also contain unique bioactive substances like D-pinitol, myo-inositol, and flavonoids, which can regulate blood glucose levels, prevent metabolic diseases, and exert antioxidant effects^[2−4]. Furthermore, M. crystallinum is a facultative crassulacean acid metabolism (CAM) halophyte, requiring partial NaCl for optimal growth^[5]. Under salt or drought stress, its photosynthetic pathway can shift from C3 to facultative CAM, achieving a water-use efficiency 3−6 times higher than typical C3 plants^[6]. Additionally, the bladder cells on its stems and leaves effectively sequester Na ions from the cytoplasm by storing natural salts, significantly mitigating ionic toxicity and osmotic imbalance caused by salt stress, further supporting its adaptive survival in high-salinity environments^[5].

Currently, global salt-affected land area exceeds 900 million hectares (accounting for approximately 10.7% of the total global land area)^[7], posing a global environmental challenge that significantly inhibits seed germination, crop development, and yield formation^[8]. To cope with salt stress, plants have evolved multi-layered regulatory mechanisms^[9]. The core of these mechanisms lies in adaptive responses achieved through the regulation of gene expression. In this process, transcription factors (TFs) play a central regulatory role at the transcriptional level by activating or repressing downstream target genes^[10].

In plants, basic helix-loop-helix (bHLH) proteins form a large TF family, second only to the MYB TF family in size^[11]. The bHLH family members are defined by the highly conserved bHLH domain, approximately 50−60 amino acids long, consisting of two functionally distinct sub-regions: an N-terminal basic region (10−15 amino acids) and a C-terminal helix-loop-helix (HLH) region (40−50 amino acids)^[12]. The basic region contains a conserved HER (His-Glu-Arg) motif that is required for binding to the E-box (CANNTG) element^[13]. Enriched in basic residues, the basic region directly recognizes and binds to DNA cis-elements like the E-box^[14] and G-box (CACGTG) target^[15]. The HLH region, composed of hydrophobic residues forming two amphipathic α-helices connected by a loop of variable length, mediates protein homo- or hetero-dimerization. This process is crucial for DNA-binding activity and functional specificity^[13].

Due to the crucial roles of bHLH genes in diverse physiological and biological processes, they have been extensively studied across various species. It has been reported that bHLH TFs play important regulatory roles in the growth and development of various reproductive tissues, including petals^[16], tapetum, and pollen^[17]. Simultaneously, some bHLH proteins participate in multiple signal transduction processes, such as light^[18] and hormone signaling^[19]. Notably, bHLH proteins are involved in responses to various abiotic stresses, such as drought^[20], salt stress^[21], low temperature^[22], and high-temperature stress^[23]. Research has shown that bHLH TFs can act as transcriptional activators in abscisic acid (ABA) signaling, regulating stomatal aperture^[19,23], while also functioning as in vivo kinase substrates for ABA, acting as stomatal opening enhancers^[24]. As a facultative CAM plant, nocturnal stomatal opening and diurnal closure are key characteristics of the CAM pathway in the M. crystallinum under salt stress. Moreover, existing studies have confirmed a close relationship between the remodeling of the photosynthetic pathway under salt stress and ABA in M. crystallinum. Therefore, in-depth research on the bHLH gene family provides an important theoretical basis for deciphering the mechanisms underlying photosynthetic pathway remodeling and salt stress tolerance in M. crystallinum.

Currently, the functions of the bHLH gene family have been extensively characterized in plants^[25], they however remain uncharacterized in M. crystallinum. Hence, this study identified 91 M. crystallinum bHLH family genes and conducted comprehensive bioinformatics analyses. The present analytical framework encompassed comprehensive characterization of the bHLH family, including family member identification, physicochemical properties, and subcellular localization prediction, phylogenetic tree construction, domain and motif analysis, gene structure analysis, gene family evolutionary tracking (gain and loss), chromosomal localization, Ka/Ks calculation, synteny analysis, co-expression network analysis, functional annotation, promoter cis-acting element analysis, expression profiling, and protein-protein interaction network analysis. The results provide a crucial reference for evolution and functional elucidation of the bHLH gene family in M. crystallinum.

Data files required for family identification of M. crystallinum and its related species Chenopodium pallidicaule, Fagopyrum tataricum, Hylocereus undatus, Chenopodium quinoa, Atriplex hortensis, Amaranthus cruentus, Simmondsia chinensis, Beta vulgaris, and Vitis vinifera were downloaded from PlantGIR (http://plantgir.cn). Simultaneously, relevant data files for Spinacia oleracea, Arabidopsis thaliana, and Oryza sativa were downloaded from Phytozome (https://phytozome.jgi.doe.gov). Identification was performed using HMMER 3.2.1 with the bHLH domain profile (Pfam, PF00010) to search for members containing the bHLH domain within the complete protein sequences of the 13 species (default parameters), using an E-value threshold of < 10⁻⁵. The domain of each protein sequence was double-verified using the CDD and PFAM databases within the InterProScan online tool (www.ebi.ac.uk/interpro/about/interproscan). Proteins containing the bHLH domain were retained as the final bHLH family members.

Based on the obtained bHLH genes, their isoelectric point (pI), molecular weight, and other molecular characteristics were predicted using the 'Protein parameter Calc (ProtParam-based)' plugin in TBtools-II software (version 2.210)^[26]. Subcellular localization of the bHLH proteins in M. crystallinum was predicted using Plant-mPLoc^[27], and WoLF PSORT^[28].

Analysis of the domain, motif, and exon-intron structure of bHLH family

The domain composition of all bHLH protein sequences was analyzed using the CDD and PFAM databases within the InterProScan online tool. Then, motifs within all bHLH protein sequences were analyzed using the MEME website (number of motifs set to 10, other parameters default)^[29]. Next, the exon-intron structure of all bHLH protein sequences was analyzed using CFVisual software (version 2.1.5)^[30]. Finally, the results were visualized using iTOL (version 7, https://itol.embl.de)^[31], and CFVisual software, respectively.

Phylogenetic tree, inference of gene gain and loss, and genotype classification of the bHLH gene family

Sequences were aligned using MAFFT-7.526^[32]. A phylogenetic tree of the bHLH family for M. crystallinum and A. thaliana was constructed using IQ-TREE2^[33]. Subsequently, a phylogenetic tree including the C3 plant M. crystallinum, the C4 plant A. cruentus, and the CAM plant H. undatus was constructed (default parameters, bootstrap number = 1,000).

To investigate the quantitative evolutionary process of the bHLH family, Notung software (version 2.9.1.5) was used to infer the history of gain and loss of bHLH family members from the ancestor to the extant species among the 13 species mentioned above, including M. crystallinum, based on previous methodologies^[34,35]. Based on previous research^[15], the bHLH family members of the study species M. crystallinum and A. thaliana were classified according to the classification, phylogenetic relationships, and topological structure established for A. thaliana.

Chromosomal localization and synteny analysis of the bHLH gene family

Chromosomal localization of the M. crystallinum bHLH gene family was performed and visualized using TBtools-II software. The Ka (nonsynonymous substitution rate) and Ks (synonymous substitution rate) values for homologous gene pairs within M. crystallinum were calculated using TBtools-II. Intra-species synteny analysis for M. crystallinum was performed using MCScanX^[36]. Additionally, two separate four-species inter-species synteny analyses were conducted: one involving M. crystallinum and the three model species A. thaliana, V. vinifera, and O. sativa; and another involving M. crystallinum and the three related species S. oleracea, F. tataricum, and H. undatus. Syntenic plots of homologous genes were generated.

Analysis of the co-expression network, functional annotation, promoter cis-acting element, and protein-protein interactions network of bHLH genes

Python scripts were used to perform Pearson correlation coefficient (PCC) analysis between the identified M. crystallinum gene set and RNA-seq expression data (FPKM values) from various sample groups (including control and salt-treated samples). Data with correlation coefficients greater than 0.95 or less than −0.95 were extracted, and co-expression network diagrams were drawn for visualization.

Functional annotation (including description, GO, nr, pfam, swissprot, and trembl) of the M. crystallinum bHLH gene family was performed using genome annotation data downloaded from TBGR (www.tbgr.org.cn). To further explore the regulatory network of bHLH genes, the 1500 bp nucleotide sequence upstream of the start codon was extracted for each M. crystallinum bHLH gene as its promoter region. Prediction of promoter cis-acting elements was then conducted using the PlantCARE website (default parameters)^[37]. Finally, the results were visualized using CFVisual software.

Furthermore, the protein-protein interaction network profile for the proteins encoded by the M. crystallinum bHLH genes was analyzed using the STRING database (Version: 12.0)^[38] (parameter setting: confidence score > 0.7).

Plant growth conditions and expression analysis

RNA-seq expression data for M. crystallinum were obtained based on sample data measured in previous literature^[6,39,40], and experimental measurements conducted by the present research group. This dataset includes FPKM expression values for all M. crystallinum genes across experimental replicates under different time points, varying salt concentration treatments, and different tissue types (root, stem, and leaf) respectively. Python scripts were then written to extract expression data for all M. crystallinum bHLH genes, and heatmaps were generated for visualization.

Ice plant seeds were imbibed at 4 °C for 2 d. They were then directly planted in soil (Klasmann-Deilmann 876) and grown at 23 ± 1 °C under 200 μmol·m⁻²·s⁻¹ white light with a 14 h/10 h day/night cycle. During the maintenance process of M. crystallinum, plants were irrigated with 0, 25, 50, 75, 100, and 125 mM NaCl every week. After 42 d, the second pair of mature leaves was collected for RNA sequencing and qRT-PCR examination.

For qRT-PCR analysis, total RNA was isolated using RNA Easy Fast Plant Tissue Kit (Tiangen, Tiangen Biotech [Beijing] Co. Ltd). First-strand cDNA was synthesized using the FastKing RT Kit (with gDNase, Tiangen, Tiangen Biotech [Beijing] Co. Ltd). qRT-PCR was performed using SuperReal PreMix Plus (SYBR Green, Tiangen, Tiangen Biotech [Beijing] Co. Ltd). One μL cDNA was used as a template. The qRT-PCR conditions were set as follows: 95 °C for 15 min, followed by 40 cycles of 95 °C at 10 s, and 60 °C at 30 s. FNR1 was used as an internal control. Primers used in the qRT-PCR are shown in Supplementary Table S1. Three replicates were employed, and the relative expression levels were calculated using the 2^−ΔΔCᴛ method. Statistical differences were determined by variance analysis. A p-value < 0.05 and an absolute fold change ≥ 2.0 signified differential expression.

A total of 91 non-redundant genes containing the typical bHLH domain were identified within the M. crystallinum genome. Comparative analysis with related species revealed 84, 82, 125, 90, 156, 153, 74, 94, and 90 bHLH genes identified in S. oleracea, C. pallidicaule, F. tataricum, A. cruentus, C. quinoa, H. undatus, B. vulgaris, S. chinensis, and A. hortensis, respectively. Additionally, 125, 115, and 158 bHLH genes were identified in the reference model species O. sativa, V. vinifera, and A. thaliana, respectively. Overall, the number of bHLH genes in M. crystallinum (91) falls within an intermediate range, being lower than quinoa, rice, V. vinifera, and A. thaliana, but higher than some related dicot species.

Further analysis showed that the 91 M. crystallinum bHLH proteins range in length from 135 to 702 amino acids. Their molecular weights range from 15.7 to 76.7 kDa, theoretical pI range from 4.47 to 9.8, instability indices range from 26.54 to 72.67, aliphatic indices range from 49.78 to 103.37, and grand average of hydropathicity (GRAVY) values range from –1.057 to –0.126, indicating a high degree of structural diversity. At the subfamily level, proteins belonging to subfamilies VII(a + b) generally possess longer amino acid sequences and larger molecular weights; conversely, subfamilies Ib(2) and IVc exhibit more similar characteristics. Notably, subfamily IIIf contains both the protein with the largest molecular weight (Mc02G01577) and the smallest molecular weight (Mc03G02933) within the entire genome, reflecting functional divergence within this subgroup (Supplementary Table S2).

Subcellular localization predictions indicated that the vast majority of M. crystallinum bHLH proteins are localized to the nucleus, consistent with their fundamental role as TFs. A small number of proteins were predicted to localize to the cell membrane, chloroplasts, or mitochondria (Supplementary Table S3). These results not only demonstrate significant diversity in the physicochemical properties and structural composition of the M. crystallinum bHLH protein family but also provide important clues for subsequent functional classification and in-depth investigation of regulatory mechanisms.

Exon-intron structure, domain, and motif analysis of the bHLH gene family

The gene structure of M. crystallinum bHLH genes were analyzed to understand their diversity. Analysis revealed that the average length of M. crystallinum bHLH genes is approximately 4,352 bp, with a maximum length of 24,507 bp, and a minimum length of 673 bp. The number of introns inserted ranges from 0 to 11. It was observed that genes such as Mc04G00185 and Mc03G02400 [subfamily III(d+e)], Mc02G01577 (subfamily IIIf), and Mc05G00505 [subfamily VII(a+b)] are relatively long, while genes in subfamilies like Ib(2) and VIIIc(2) are relatively short. Untranslated regions (UTRs) flank the gene sequence in most subfamilies; however, subfamilies such as VIIIb, VIIIc(1), and VIIIc(2) lack UTRs. Subfamilies like Ia and Ib(2) are intron-poor (≤ three introns per gene), whereas subfamilies IVb, IVc, and IVd are intron-rich (> three introns per gene). Furthermore, genes within subfamilies III(d+e), VIIIb, and Orphans lack introns entirely. This indicates functional diversity and complexity among these genes (Fig. 1, Supplementary Table S4).

Figure 1. 

Motif and gene structure analysis of the M. crystallinum bHLH gene family. (a) The phylogenetic tree of 91 M. crystallinum bHLH. (b) Conserved motif distribution of bHLH protein in M. crystallinum. (c) The genetic structure of the bHLH gene in M. crystallinum. The information of each gene structure is provided in Supplementary Table S4. The red rectangle and the yellow rectangle represent the untranslated region (UTR), and the exon respectively. Introns are indicated by black lines. The scale is shown at the bottom of the chart.

The domain composition of the M. crystallinum bHLH gene family was also analyzed, revealing diversified domain structures. The results showed 16 distinct types of bHLH domains with high diversity. Protein members within the same subfamily generally contain similar types of bHLH domains, with a few exceptions [e.g., Ib(2), IIIf] (Supplementary Fig. S1).

Regarding motif composition, 10 motifs were identified among the core M. crystallinum bHLH proteins. Motif 1 and Motif 2 are ubiquitous across all sequences. Subfamily XII contains Motif 5; additionally, protein Mc06G01080 [subfamily VII(a+b)], Mc04G00224 [subfamily III(d+e)], and Mc01G00955 (subfamily Vb) also contain Motif 5, while others do not. Most proteins in subfamilies III(d+e) and Orphans contain Motif 6 and Motif 9, which can also be found in a few other subfamilies, such as protein Mc06G00350 in subfamily XIV. Overall, while motif composition varies, proteins belonging to the same evolutionary clade tend to be similar. Almost all M. crystallinum bHLH proteins lack some motifs, and duplication events are relatively infrequent (Fig. 1).

Phylogenetic tree and gain/loss in the bHLH gene family

To achieve a detailed classification of the M. crystallinum bHLH gene family, a phylogenetic tree including bHLH family genes from M. crystallinum and the model species A. thaliana were constructed (Fig. 2). The maximum likelihood (ML) phylogenetic tree, constructed based on 91 M. crystallinum bHLH protein sequences and A. thaliana reference sequences, revealed that the family can be classified into 27 subfamilies: Ia, Ib(1-2), II, III(a+c), IIIb, III(d+e), IIIf, IVa-d, Va, Vb, VII(a+b), VIIIa-b, VIIIc(1), VIIIc(2), IX, X, XI, XII, XIII, XIV, XV, and Orphans. This classification was based on topological structure and the evolutionary relationships of homologous genes in A. thaliana (Fig. 2). The phylogenetic tree showed that most minor subfamilies belonging to the same major subfamily clustered together on the same branch and were closely related. Furthermore, M. crystallinum genes belonging to subfamilies such as Ib(2), XII, and III(d+e) were relatively more numerous.

Figure 2. 

Phylogenetic analysis of the bHLH gene family from M. crystallinum and A. thaliana. Phylogenetic trees were constructed using the bHLH protein sequences and different colors indicate different clades.

By reconciling the phylogenetic tree with the species tree, Notung software was used to infer the history of gene gain and loss events in the bHLH gene family across the 13 study species, including M. crystallinum (Fig. 3). The results indicated that the current number of bHLH genes in M. crystallinum (91) was formed from its ancestral genes through the loss of 66 genes and the duplication of seven genes. The evolutionary paths for other species were as follows. The common ancestor of the 13 study species contained 135 bHLH genes. It experienced the gain of 57 genes and the loss of 67 genes, forming the extant bHLH family in the monocot model species O. sativa. This ancestor experienced the gain of 68 genes and the loss of 12 genes, resulting in an ancestral pool of 191 genes for the remaining 12 species. This ancestral pool then underwent the gain of 64 genes and the loss of 97 genes, forming the extant bHLH family size in the dicot model species A. thaliana. The ancestral pool of 191 genes for the 12 species experienced the gain of 21 genes and the loss of 28 genes, forming an ancestral pool of 184 genes for V. vinifera, and the other 10 species. This pool subsequently experienced the gain of 24 genes and the loss of 93 genes, forming the extant bHLH family size in the genomic model species V. vinifera. Similarly, the extant bHLH family sizes in the related species of M. crystallinum were formed through gene loss and duplication events.

Figure 3. 

Analysis of the gain and loss of the bHLH gene and evolutionary time in 13 species. Schematic diagram of the gain and loss of the bHLH gene in M. crystallinum and their closely related species. The green rectangle and the numbers in yellow font on the right side of the species picture represent the quantity of bHLH genes in ancestors and existing species. The + sign and the - sign respectively indicate the gain and loss of genes.

Through comparison, it was found that A. thaliana experienced the highest number of bHLH gene losses (97), while the related species C. quinoa experienced the fewest losses (10). Notably, among these 13 species, only C. quinoa had more gene duplications (23) than losses (10). In all other species, gene duplication was less frequent than gene loss. Overall, during the evolution of the bHLH family in the 13 study species, gene loss (1,046 events) significantly outnumbered gene duplication (733 events) (Fig. 3).

Synteny analysis of bHLH genes in M. crystallinum and related species

To investigate the chromosomal distribution characteristics of the M. crystallinum bHLH gene family, chromosomal localization of all 91 bHLH genes in M. crystallinum and their orthologs in three closely related species were performed (Fig. 4). The results showed that M. crystallinum bHLH genes were distributed across all nine chromosomes, with significant non-uniform distribution patterns observed in both M. crystallinum and the three related species. In M. crystallinum, chromosomes 2, 3, and 4 exhibited dense gene clustering, while chromosome 7 contained only a single gene (subfamily IVc). However, none of the three related species showed such extreme cases of chromosomes containing very few genes. Overall, most M. crystallinum bHLH genes were located near chromosome termini, a distribution pattern consistent with their orthologs in related species.

Figure 4. 

Chromosome localization of the bHLH gene family in M. crystallinum, A. thaliana, H. undatus, and S. oleracea. The left scale indicates the chromosome length (Mb), with bHLH gene markers on the right side of each chromosome. Among them, blue represents M. crystallinum, green represents A. thaliana, purple represents H. undatus, and orange represents S. oleracea.

Expansion of the M. crystallinum bHLH gene family primarily occurred through dispersed duplication and segmental duplication (Supplementary Table S5). Intraspecific synteny analysis using MCScanX identified eight synteny gene pairs distributed across chromosomes 3–9 (Fig. 5). Exon-intron structures of most synteny pairs were highly conserved (only Mc03G00602-Mc03G02810 differed in intron category), indicating functional conservation post-duplication.

Interspecific synteny analysis revealed evolutionary relationships between M. crystallinum and other species. Synteny gene pairs with A. thaliana and O. sativa were predominantly concentrated in terminal chromosomal regions (Fig. 5c), whereas those with the closely related species H. undatus covered broader segments, implying shared ancestral genomic architecture. Ka/Ks analysis showed all homologous gene pairs had Ka/Ks ratios < 1 (Table 1), indicating purifying selection during evolution.

Co-expression, protein network, and cis-acting element of M. crystallinum bHLH genes

Co-expression analysis of transcriptome data extracted gene pairs with PCC > 0.95, revealing exclusively positive correlations (all PCC > 0.95; none < −0.95; Fig. 6). Network analysis demonstrated that bHLH genes correlated highly with diverse genes, highlighting their central transcriptional roles and synergy with other families. Mc03G00604 co-expressed significantly with MYB family gene Mc01G01299 (MYB20, lignin metabolism), glycosyltransferase gene Mc01G00687 (mucilage biosynthesis), 2OG-Fe(II) oxygenase gene Mc09G00173 (disease defense), and bHLH gene Mc03G00603 (bHLH126). Mc05G01224 co-expressed with WRKY family gene Mc04G02119 (WRKY65).

Figure 6. 

Co-expression network of the bHLH gene family in M. crystallinum. Based on the correlation of expression levels among the genes of M. crystallinum, a co-expression network is formed. The light purple circle indicates that the genes M. crystallinum do not belong to the bHLH family, and the genes belonging to the bHLH family are marked with different colors. The correlations in the figure are all greater than 0.95.

In the intraspecific synteny pair Mc03G00604-Mc03G00603 (PCC = 0.992), Mc03G00603 was annotated as transcription factor (TF) bHLH120 (AT5G51790), with GO enrichment for 'DNA-binding TF activity (GO:0003700)' and 'RNA polymerase II transcription regulation (GO:0006357)', suggesting a core role in transcriptional cascades. Co-expression of bHLH genes with ABC transporter (Mc02G02356), GST (Mc03G01764), and copper transporter (Mc03G00444) genes further implicated roles in transport and detoxification.

Notably, bHLH members showed significant co-expression with key ABA/JA pathway genes: Mc01G00468 and Mc02G01019 co-expressed with the ABA signaling gene Mc06G01384. Mc03G00601 co-expressed with the ABA receptor gene Mc08G01193. Mc03G00601, Mc03G00603, and Mc03G00604 co-expressed with the JA metabolic gene Mc02G00882. Mc03G00600, Mc03G00603, and Mc08G01910 co-expressed with the JA-responsive TF gene Mc06G02065. Mc03G00601 co-expressed with both ABA receptor and JA metabolic enzyme genes, indicating potential cross-talk mediation. Additionally, Mc03G01764 is strongly associated (PCC = 0.989) with glutathione S-transferase GSTU8 (detoxification). Genes linked to Mc02G01019 (Mc05G01832, Mc06G01713) showed synergy with auxin transporter PIN2 and ethylene response factor WIN1 (Supplementary Table S6).

Annotation analysis of the M. crystallinum bHLH TF family identified key regulatory members: Mc09G00905 (PIF4) involved in light signaling, Mc02G01681 (PIL5/PIF1) associated with circadian regulation and cold response, and Mc04G00224 (SCRM2/ICE2) functioning in negative gene expression regulation. GO analysis revealed enriched molecular functions including DNA binding (GO:0003677), TF activity (GO:0003700), and protein dimerization (GO:0046983), with significant involvement in biological processes such as light response (GO:0009416), seed germination regulation (GO:0010114), and iron ion homeostasis (GO:0055072). Mc06G01874 was annotated as a core regulator of cellular iron starvation response (GO:0010106). Domain analysis confirmed the universal presence of the conserved HLH domain (PF00010.27) in all members, with additional bHLH-MYC_N domains (PF14215.7) detected in subsets (e.g., Mc02G01577), implying MYC-like regulatory networks. Functional complexity was noted: Mc03G02400 (MYC2-like) contains both HLH domain and metabolic enzyme annotation (2-alkenal reductase), implying dual regulatory-metabolic roles (Supplementary Table S7).

Analysis of 1.5 kb promoter sequences upstream of start codons identified enriched cis-acting elements (Fig. 7): (1) Stress-responsive elements: LTR (low-temperature), MBS (drought), TC-rich repeats (defense), with salt-specific elements ARE and GC-motif enriched in Mc05G01869 and Mc06G01874; (2) Hormone-responsive elements: ABRE (ABA), CGTCA/TGACG-motif (MeJA), GARE/P-box (GA); (3) Light-responsive elements: G-box (CACGTG), I-box; (4) Other functional elements: circadian (diurnal rhythm), MSA-like (cell cycle). Stress-responsive elements predominated, with notable clustering (e.g., LTR in Mc03G03147) (Fig. 7). These results indicate M. crystallinum bHLH genes integrate light, hormone, and stress signals through complex promoter cis-element networks (Supplementary Tables S8–S10).

Figure 7. 

Statistics of cis-acting elements in the bHLH promoter of M. crystallinum. Red represents light response, blue represents stress response, green represents hormone response, and purple represents others.

To understand the biological functions of bHLH proteins, the STRING database was used to analyze the protein interaction networks of their expressed proteins. The results showed that 89 candidate proteins of M. crystallinum successfully match the homologous proteins of A. thaliana, with sequence similarities ranging from 25.9% to 91.1%. The core functional modules include: (1) MYC2 can bind to cis-elements such as G-box and Z-box, and participate in the signaling pathways of light, ABA, and jasmonic acid (JA). PIF1 activates genes such as SOM by binding to the G-box cis-element, and cooperates with bHLH122 to enhance the functions mediated by JA/ABA signals. (2) The stomatal dynamic regulation module: SPCH and FAMA interact with SCRM to inhibit excessive stomatal closure. BHLH85 and GL3 are involved in the life processes of epidermal cells and negatively regulate stomatal formation. (3) The metabolic reprogramming module: BIM1 interacts with PIF4 and PIF3 to activate photoperiod response genes. TT8 and BHLH93 enhance cell membrane stability by regulating flavonoid synthesis genes. Members such as BHLH160 form a specific interaction network with BHLH162 in M. crystallinum, possibly optimizing the expression of osmotic adjustment genes through variations in the bHLH domain (Fig. 8). This network reveals that the bHLH family supports the adaptive evolution of M. crystallinum in high-salt environments through both conserved stress signal transduction and species-specific metabolic integration.

Figure 8. 

Protein interaction network of the bHLH gene family in M. crystallinum. PPI network based on the bHLH homolog in A. thaliana.

Expression analysis profiling of bHLH under salt stress in M. crystallinum

Based on the expression data analysis of M. crystallinum, the gene expression of M. crystallinum shows significant tissue preference (Supplementary Table S11): Mc08G00451 is highly expressed in leaves (40) and stems (47), but is relatively low in roots (6), while Mc03G02400 is highly expressed in roots (138) and stems (105), but is relatively low in leaves (47). Mc09G00237 maintains high expression in roots (26.72), stems (22.67), and leaves (18.27), possibly participating in the multi-tissue coordinated regulation of basic metabolic pathways. Mc04G00255 has a significantly higher expression in leaves (18) compared to roots (2), and stems (0) (Fig. 9a), suggesting that it may be preferentially involved in stress signal transduction in leaves. The expression of Mc06G00859 gene increases first, then slightly decreases, and then stabilizes as the salt concentration continuously rises, indicating that it may be expressed under salt conditions.

Figure 9. 

Expression pattern analysis of bHLH genes in M. crystallinum. Expression levels of the bHLH gene in the (a) roots, stems and leaves. (b) M. crystallinum grown under different salt concentrations (0, 25, 50, 75, 100, and 125 mM NaCl). (c) M. crystallinum with control and salt treatment (500 mM NaCl) collected at 8:00 and 20:00 from 5 to 7 d after treatment. (d) M. crystallinum under different salt treatments (0, 140, 250, and 500 mM NaCl). (e) M. crystallinum treated with 500 mM NaCl salt concentration and the control group on the first and second days. (f) M. crystallinum treated with CK, LS (100 mM), MS (300 mM), and HS (500 mM) NaCl. The expression levels increase from blue to white and then to red.

The expression of multiple genes in M. crystallinum shows significant dynamic responses with the increase of salt concentration. Mc05G01869 increases to 33 under T75mM treatment, approximately 7.5 times higher than the control group (4.38), suggesting that it participates in the response under high salt concentration. Mc03G00601 has an expression level of 0.13 in the control group, but significantly decreases to 0 under salt treatment (T125mM), suggesting that its function may be inhibited by salt stress (Fig. 9b).

Furthermore, the expression levels of M. crystallinum genes vary in different experiments. For example, Mc06G01846 shows a trend of increasing first and then decreasing under salt treatment (500 mM NaCl) (Fig. 9c), but in the data of four different concentrations (0, 140, 250, 500 mM), the expression level increases with the increase of salt concentration (Fig. 9d). When M. crystallinum is grown for 42 d, genes such as Mc04G00301, Mc05G01869, Mc04G02348, and Mc03G02400 show a decreasing trend on the second day compared to the first day under 500 mM NaCl salt concentration treatment (Fig. 9e). The expression levels of M. crystallinum genes under different concentrations (CK, LS: 100 mM NaCl, MS: 300 mM NaCl, HS: 500 mM NaCl) show a continuous increase (Mc06G02041), first increase and then decrease (Mc05G01869, Mc09G01047), and decrease with the increase of concentration (Mc05G00898), presenting different expression trends (Fig. 9f).

To validate the response of M. crystallinum bHLH transcription factors to salt stress, eight significantly differentially expressed bHLH genes were selected from transcriptome data for q-PCR verification; the results revealed that Mc05G01869 [subfamily III(a+c)], Mc03G02400 and Mc04G00185 [subfamily III(d+e)], and Mc06G00859 (subfamily Vb) all exhibited an initial increase followed by a decrease in expression levels as salt concentration increased, while Mc08G01331 (subfamily Vb) demonstrated a gradually increasing trend in relative expression (Fig. 10).

Figure 10. 

qRT-PCR analysis of eight differential expression bHLH gene in leaf under different concentrations of salt irrigation. Different letter designations for the same gene indicate significant differences (p < 0.05).

Mounting evidence indicates that plant bHLH genes participate in diverse physiological and biological processes, including development and stress responses^[41]. However, research on the facultative CAM plant M. crystallinum, renowned for its extreme salt tolerance, remains limited. The release of a high-quality M. crystallinum genome^[6] enables in-depth investigation of its bHLH family. Based on A. thaliana phylogenetic classification, bHLH genes from M. crystallinum, A. thaliana, and closely related species, were categorized into 27 subfamilies, with M. crystallinum distributed across 23 subfamilies. Thus, a comprehensive classification and systematic analysis of the M. crystallinum bHLH family was conducted.

Evolution of the M. crystallinum bHLH gene family

bHLH genes emerged in land plants 440 million years ago^[15] and expanded dramatically in higher plants^[42]. In this study, expansion of the M. crystallinum bHLH family primarily occurred through dispersed duplication (56 genes), followed by whole-genome/segmental duplication (16 genes), tandem duplication (TD) (15 genes), and singleton (four genes). This indicates dispersed duplication as the dominant expansion mechanism, consistent with findings by Yu et al.^[43]. Polyploidization events significantly drive family expansion^[44], which was validated in M. crystallinum. Notably, TD contributed to bHLH family expansion at a level comparable to genome doubling. Previous studies demonstrated that TDs facilitated functional innovation in disease resistance genes^[45], C4 photosynthesis genes^[44], and aspartate protease genes^[46]. The present findings revealed that M. crystallinum bHLH members evolved under multiple duplication mechanisms, providing insights for future studies on evolution-function relationships.

The Ka/Ks ratios of all synteny duplicated gene pairs in M. crystallinum were less than 1 (e.g., Mc04G02580–Mc07G01674 pair, Ka/Ks = 0.11), indicating that the bHLH genes in M. crystallinum have undergone strong purifying selection during evolution (Table 1). By exploring the synteny of bHLH genes with other species, a significant number of synteny results were identified, suggesting the evolutionary conservation of the bHLH gene family in M. crystallinum and other higher plants. M. crystallinum formed 61, 50, and 120 synteny gene pairs with its close relatives S. oleracea, F. tataricum, and A. cruentus, respectively, with the highest number of synteny gene pairs observed with A. cruentus, implying a closer evolutionary relationship between the M. crystallinum bHLH family and A. cruentus.

Structural analysis of the M. crystallinum bHLH family indicates that genes within the same phylogenetic tree classification exhibited overall structural similarity, consistent with findings from most gene family studies^[43]. Additionally, previous studies have demonstrated that AtbHLH106 (At2g41130) responds to salt stress^[47], its ortholog in M. crystallinum is Mc06G01846, which clusters within the same phylogenetic clade, and exhibits significant upregulation under salt stress. However, divergent responses were observed for certain M. crystallinum genes (Mc01G00860, Mc08G00108, Mc05G00898) phylogenetically close to the salt-responsive A. thaliana genes AtbHLH121 (At3g19860) and AtbHLH122 (At1g51140), as these ice plant homologs showed minimal salt induction. Furthermore, salt-tolerant A. thaliana genes such as AtbHLH092 (At5g43650) and AtbHLH112 (At1g61660) are evolutionarily distant from M. crystallinum bHLH members, reflecting functional divergence between species.

The number of introns in M. crystallinum bHLH genes varied significantly (0–11), with subfamilies III(d + e) and VIIIb being intronless, while subfamilies IVb, IVc, and XII were intron-rich (Supplemental Table S4)^[48]. This structural divergence may enhance functional diversity through alternative splicing mechanisms. Focused on the structural differences in synteny gene pairs, it was found that most synteny gene pairs within M. crystallinum belonged to the same structural category (e.g., synteny gene pairs Mc03G00600–Mc03G02811 belonged to the intron-poor category, Mc05G00898–Mc06G02047 belonged to the intron-rich category, while Mc03G03147–Mc04G00224 had no intron insertions). However, there was one exception: the synteny gene pair Mc03G00602–Mc03G02810, where the former was intron-poor (two introns) and the latter was intron-rich (nine introns).

More deeply, the evolutionary trajectory of bHLH family size changes in M. crystallinum and its close relatives was inferred. In the context of the bHLH family across 13 examined species, it was deduced that the ancestral bHLH family size was at least 135 genes. Overall, the number of bHLH gene losses far exceeds the number of gene duplications. Previous studies have shown that whole-genome duplication (WGD) events are often accompanied by substantial gene loss^[49], and the bHLH family is no exception. In fact, similar phenomena have been observed in studies of other gene families. Intriguingly, during the duplication and loss events at ancestral nodes, some nodes showed a higher number of duplications than losses. For example, the ancestral bHLH family of M. crystallinum, H. undatus, A. cruentus, B. vulgaris, S. oleracea, A. hortensis, C. quinoa, and C. pallidicaule experienced more duplications than losses (46:30). This reflected the distinct evolutionary histories of bHLH families in various plants. Therefore, the present study elucidated the evolutionary history of the bHLH family in M. crystallinum and comparative species, providing a valuable reference for studying gene duplication, loss, and retention in other species.

Function of bHLH in M. crystallinum and exploration under salt stress

Expression profiling was conducted and validated through qRT-PCR experiments showing that specific M. crystallinum genes exhibit significant differential expression under salt stress. The involvement of bHLH TFs in M. crystallinum salt stress response aligns with the findings of Takahashi et al.^[24], who reported that the bHLH TF Mcr010456.002 (At2g42280.1) in M. crystallinum participates in stomatal movement by activating genes encoding inward-rectifying potassium channels. Guan et al.^[50] also found that M. crystallinum bHLH TFs may be involved in mechanisms of osmotic regulation and stomatal regulation. It is noteworthy that Mc03G02400, a MYC2-like transcription factor, likely increases salt sensitivity by regulating proline biosynthesis and the JA signaling pathway^[51]. PIF (Phytochrome-Interacting Factors), belonging to the bHLH transcription factor superfamily, serves as a central hub for light signal transduction and modulating ABA signaling pathways^[52]. Mc09G00905 is involved in light signaling, and Mc02G01681 is associated with circadian regulation, which may be related to the function of regulating photosynthetic pathways and enhancing salt tolerance.

Co-expression network analysis provided corroborating evidence: Mc05G01224 showed significant co-expression (PCC = 0.989) with glutathione transferase GSTU8 (Mc03G01764), and Mc04G00225 displayed strong co-expression (PCC = 0.999) with calcium exchanger CAX9 (Mc04G01942) (Fig. 8), suggesting that bHLH genes may enhance plant stress tolerance by coordinating ion homeostasis (e.g., Na⁺ efflux) and antioxidant defense^[53]. Additionally, Mc03G02810, the longest bHLH gene in M. crystallinum (Supplementary Table S4), exhibits high correlation with multiple genes. Among those co-expression genes, Mc01G00538 and Mc04G01204 are involved not only in gene silencing regulation but also in responses to salt stress, suggesting related functions for Mc03G02810. The qRT-PCR verification of the bHLH gene (Fig. 10) also confirmed our correlation analysis: the expression dynamics of Mc06G01846 showed upregulation when salinity increased, which was consistent with the effect of A. thaliana genes. In summary, the present work reveals the response pattern of ice plant bHLHs under salt stress. However, more in-depth research is needed to elucidate its specific salt-stress response mechanisms.

This study systematically analyzed the molecular characteristics and functional mechanisms of the bHLH gene family in the facultative CAM plant M. crystallinum through genome-wide analysis. A total of 91 bHLH genes were identified in its genome, revealing that gene loss dominates the family's evolution, with dispersed duplication as the primary driving force for expansion, and synteny gene pairs within the species generally undergoing purifying selection. The promoter regions are enriched with stress-responsive elements and hormone-responsive elements, suggesting that bHLH genes regulate salt adaptation by integrating light signals, ABA/JA pathways, and redox balance. Co-expression network and protein interaction analyses revealed that M. crystallinum bHLH genes collaborate with MYB, GST, and other gene families to form a multidimensional regulatory network, with gene expression under salt stress exhibiting tissue specificity. This study provides important theoretical insights into elucidating the functions and salt tolerance mechanisms of the M. crystallinum bHLH gene family.