Genome-wide identification and comprehensive characterization of the ABF/AREB transcription factor family in Medicago sativa L. and their responsiveness to abiotic stress

Alfalfa (Medicago sativa L.) is one of the most important forage legumes worldwide, serving as a vital source of protein in animal feed^[1] and playing a significant role in soil improvement and sustainable agricultural ecosystems^[2]. However, abiotic stresses such as drought, salinity, and cold frequently limit alfalfa's growth and productivity^[3]. Therefore, elucidating the molecular mechanisms underlying its stress responses and identifying key regulatory genes are essential for the genetic improvement of stress tolerance in alfalfa^[4].

In response to environmental challenges, plants activate complex signaling cascades that regulate the expression of downstream functional genes, thereby enhancing their adaptability^[5]. Abscisic acid (ABA) is a central phytohormone that mediates plants' adaptation to abiotic stresses such as drought, salinity, and cold^[6]. Within the ABA signaling pathway, the ABF (ABA-responsive element binding factor) family, also known as AREB (ABA-responsive element binding protein), represents a key group of transcription factors that serve as critical regulators in ABA-dependent stress responses^[7]. ABF proteins specifically recognize and bind to ABA-responsive elements (ABREs) in the promoters of stress-related genes, activating their expression and ultimately improving plants' stress tolerance^[8].

With advances in plant genomics and bioinformatics, the ABF gene family has been systematically identified and functionally characterized in multiple plant species^[9]. For instance, in Arabidopsis thaliana, ABF members such as AtABF1, AtABF2, AtABF3, and AtABF4 have been shown to participate in drought and salt stress responses, with their overexpression significantly enhancing stress tolerance^[10−12]. Similarly, in rice (Oryza sativa), genes such as OsABF1 and OsABI5 play important roles in ABA signaling and stress adaptation^[13]. Collectively, these studies indicate that the ABF family is evolutionarily conserved and functionally crucial in plants' stress responses.

Although the ABF family has been extensively studied in model plants (e.g., Arabidopsis) and some crops, a dedicated genome-wide analysis of this family in alfalfa – an agronomically important forage species – remains lacking. As an autotetraploid species, alfalfa possesses a complex genome structure that potentially harbors a larger and more functionally diverse ABF gene family than diploid plants^[14]. The release of high-quality genome sequences, such as the 'Xinjiangdaye' assembly^[15], provides a valuable foundation for the systematic identification and characterization of gene families in alfalfa at the whole-genome level.

In this study, we aimed to perform a genome-wide identification of the ABF gene family in alfalfa, and analyze their chromosomal distribution, gene structure, conserved domains, phylogenetic relationships, and cis-regulatory elements. Furthermore, by integrating publicly available transcriptomic datasets^[16,17] with quantitative real-time polymerase chain reaction (qRT-PCR) validation, we examined the expression patterns of MsABF genes across different tissues and under various abiotic stresses, including drought, salinity, and cold. This study is expected to provide insights into the regulatory roles of the ABF family in ABA-mediated stress responses in alfalfa and to identify potential candidate genes for molecular breeding aimed at enhancing stress tolerance.

Members of the ABF/AREB gene family in M. sativa were identified through an integrated homology and domain-based search strategy. First, the protein sequences of nine known Arabidopsis ABF/AREB members (AtABF1, AtABF2/AREB1, AtABF3, AtABF4/AREB2, AtABI5, AtDPBF2, AtAREB3, AtDPBF4, and AtbZIP15) were retrieved from The Arabidopsis Information Resource (TAIR, www.arabidopsis.org) and used as queries to perform a BLASTP search against the M. sativa proteome (obtained from https://figshare.com) with a stringent e-value cutoff of 1e⁻⁵. Concurrently, the hidden Markov Model (HMM) profile for the conserved bZIP domain (PF00170) was downloaded from the Pfam database (http://pfam.xfam.org) and used to scan the same proteome using HMMER 3.0 (http://hmmer.org) with the default parameters. Candidate sequences obtained from both the BLASTP and HMMER searches were pooled, and redundant sequences were removed. The presence of the definitive bZIP domain in each unique candidate was further verified using the NCBI Conserved Domain Database (CDD, www.ncbi.nlm.nih.gov/Structure/cdd). Proteins that passed this validation were designated as MsABF family members. The same pipeline was applied to identify ABF/AREB family members in Glycine max (soybean). The basic physicochemical properties of the confirmed MsABF proteins, including molecular weight, theoretical isoelectric point (pI), instability index, and other relevant parameters were calculated from their amino acid sequences using TBtools software.

Phylogenetic and structural analysis of the ABF/AREB family

To investigate the evolutionary relationships and structural characteristics, a phylogenetic tree was constructed using the full-length protein sequences of MsABF members and their Arabidopsis orthologs. Multiple sequence alignment and tree inference were performed in MEGA (version 11) using the maximum likelihood method with 1,000 bootstrap replicates. Conserved protein motifs were identified via the MEME suite with the number of motifs set to 10. Gene structure information (exon/intron organization) was extracted from the genome annotation file. All results, including the phylogenetic tree, motif distribution, and gene structure, were integrated and visualized using TBtools-II.

Chromosomal localization, gene duplication, and collinearity analysis

The chromosomal positions of the identified MsABF genes were extracted from the M. sativa genome annotation (GFF3 file) and visualized. To investigate gene duplication events, tandem duplications were defined as homologous gene pairs located within 200 kb on the same chromosome. Segmental duplications and collinearity relationships were analyzed using MCScanX, with a BLASTP E-value cutoff of 1e⁻⁵. The analysis included both intragenomic collinearity within M. sativa and intergenomic comparisons with two reference species, Arabidopsis and Glycine max, to elucidate the evolutionary patterns. All results were integrated and visualized using TBtools.

Analysis of oromoter cis-acting elements

To investigate the potential regulatory mechanisms of MsABF genes, the 2.0 kb genomic sequences upstream of the transcription start site for each MsABF gene were extracted as putative promoter regions. These sequences were submitted to the PlantCARE database for in silico identification of cis-acting regulatory elements. The predicted elements were categorized and visualized using TBtools.

Expression profiling analysis using public RNA-seq data

To examine the expression patterns of MsABF genes, publicly available RNA-seq data were retrieved from the NCBI Sequence Read Archive (SRA) database. These datasets encompassed six distinct tissues (roots, elongated stems, pre-elongated stems, leaves, flowers, and nodules; BioProject: SRP055547) and samples subjected to salt, drought, and cold stress treatments (SRA accessions SRR7091780–SRR7091794 and SRR7160313–SRR7160357). The expression levels of MsABF genes were extracted and visualized using TBtools to assess the specificity of their spatial expression and their stress responsiveness.

Prediction of the protein–protein interaction network

A protein–protein interaction (PPI) network for the MsABF proteins was predicted de novo using the STRING database (version 11.5), with M. sativa selected as the target organism and a minimum required interaction score > 0.700. The raw network data obtained from STRING was subsequently imported into Cytoscape (version 3.9.1) software for network construction, visualization, and topological analysis.

Plant materials, stress treatments, and RT-qPCR analysis

Seeds of the alfalfa cultivar 'Zhongmu No. 4', obtained from the Institute of Animal Science, Chinese Academy of Agricultural Sciences, were germinated on moist filter paper. Seedlings with radicles reaching 4–5 cm were transferred to a nutrient solution and grown under controlled conditions (16/8 h light/dark cycle, 24/20 °C day/night temperature, 70%–80% relative humidity) for 21 d. Uniform 21-d-old plants were then subjected to abiotic stress treatments: Salt stress (200 mM NaCl), osmotic stress (15% polyethylene glycol [PEG]-6000), and cold stress (4 °C). Untreated plants were maintained as controls. Leaf samples (the third fully expanded trifoliate leaf from the apex) were collected at 0, 6, 12, and 24 h after the onset of each treatment, with each time point having three biological replicates.

Total RNA was extracted from all samples and reverse transcribed to obtain cDNA. Gene-specific primers for MsABF genes were designed using Primer 5.0 software. The RT-qPCR assay was performed on an Applied Biosystems 7,500 system using SYBR Premix Ex Taq, with three technical replicates per biological sample. The MsActin (EU664318) gene served as the internal reference for normalization, and relative gene expression levels were calculated using the 2^−ΔΔCᴛ method.

A comprehensive genome-wide search identified a total of 46 MsABF genes in the alfalfa genome, each containing the characteristic bZIP domain (PF00170). These genes were systematically designated MsABF1 to MsABF46 according to their sequential order along the chromosomes from chromosome (chr) 1.1 through to chr 8.3, arranged in a left-to-right manner (Fig. 1). Chromosomal localization revealed an uneven distribution across the 27 chromosomes (chr 1.1 to chr 8.3), with the highest density observed on chr 4.4, which harbors five members (MsABF21 to MsABF25).

Figure 1. 

Chromosome distribution of ABF genes in M. sativa.

The gene IDs and corresponding protein sequences of all 46 MsABF members are listed in Supplementary Table S1. The encoded proteins exhibited diverse physicochemical characteristics, as detailed in Supplementary Table S2. Specifically, the protein lengths ranged from 156 to 688 amino acids, with molecular weights varying between 17,336.19 and 74,994.38 Da. Theoretical isoelectric points (pI) ranged from 4.92 to 9.74, whereas instability indices predicted that the MsABF17 members were stable. All MsABF proteins showed negative grand average of hydropathy (GRAVY) values, confirming their hydrophilic nature.

Phylogenetic and structural analysis of the MsABF family

To elucidate the evolutionary relationships among the MsABF proteins, a phylogenetic tree was constructed with their Arabidopsis orthologs. The analysis resolved the MsABF family into seven distinct groups (Groups I−VII) (Fig. 2). Group I was the largest, containing 15 members, whereas Groups III and V each comprised 9 members. Groups II, IV, and VI contained intermediate numbers, and Group VII was the smallest with only two members.

Figure 2. 

Phylogenetic tree of ABF proteins in M. sativa and Arabidopsis. The blue dots represent M. sativa, and the pink dots represent Arabidopsis.

An analysis of the conserved protein motifs identified a total of 10 distinct motifs (Fig. 3a, detailed consensus sequences in Fig. 3d). Motif 2 was present in all members, likely corresponding to core functions such as DNA binding. In contrast, Motif 8 was present only in a subset of members (e.g., MsABF11, MsABF13), suggesting its involvement in specialized regulation. The extreme cases of MsABF29 (containing only a single motif) and MsABF9 (containing all 10 motifs) imply a trend toward subfunctionalization within the family. Domain analysis confirmed that all members possess the canonical bZIP domain (Fig. 3c). Gene structure was highly conserved, with all members lacking UTRs (untranslated regions) and exhibiting relatively conserved intron numbers (Fig. 3b). This structural consistency suggests that functional diversification within the MsABF family in alfalfa likely occurred primarily through gene duplication rather than structural rearrangement.

Figure 3. 

(a) The motif structure diagram of ABF genes. (b), (c) The basic structure of ABF genes. (d) Basic information of ABF protein motifs in M. sativa.

Gene duplication and collinearity analysis

Genomic analysis revealed distinct expansion patterns within the MsABF family. As shown in Fig. 4a, 53 segmental duplication events in total were identified, with chr 3.1 and chr 3.4 showing the highest duplication density (8 events each). Notably, three genes, MsABF2, MsABF18, and MsABF20, were each involved in six duplication events, suggesting their particular importance in recent evolutionary expansion and potential functional diversification.

Figure 4. 

(a) Schematic diagram of the linear distribution of ABF genes in M. sativa. Segmental replication events are marked in red. The red and blue stripes in the inner ring indicate the genes' density on the chromosome. (b) Collinearity of MsABF genes with Arabidopsis and Glycine max genes. The blue line represents homology between the genes.

Interspecific collinearity analysis revealed extensive synteny between M. sativa and related species. In total, 28 collinear gene pairs between M. sativa and A. thaliana, and 86 pairs between M. sativa and G. max (Fig. 4b) were found. These correspond to 16 and 36 MsABF genes exhibiting synteny with A. thaliana and G. max, respectively. Notably, nine MsABF genes displayed no detectable collinearity with either reference species, suggesting that they may represent lineage-specific genes or that they have undergone substantial genomic rearrangement.

Analysis of promoter cis-acting elements

Promoter analysis of the MsABF gene family revealed complex regulatory architectures potentially governing their expression patterns (Fig. 5b−d). In total, 826 major cis-acting elements were identified in the promoter regions 2.0 kb upstream. These elements were categorized into four primary functional groups: Growth and development-related (2 types, 10 elements), light-responsive (6 types, 346 elements), MYB-binding sites (2 types, 60 elements), stress-responsive (5 types, 184 elements), and hormone-responsive elements (5 types, 226 elements). The distribution of these elements across different MsABF members showed significant phylogenetic correlations (Fig. 5c). Light-responsive elements were predominantly concentrated in members from phylogenetic Groups I and III. Stress-responsive elements (drought, low-temperature, and wounding response elements) showed higher density in Groups I and VI. Similarly, hormone-responsive elements (ABA, jasmonic acid [JA], and salicylic acid [SA]-responsive elements) were primarily enriched in Groups I and III. This cis-element landscape reveals distinct regulatory specializations among MsABF subfamilies. Members of Groups I and III appear to integrate multiple environmental signals, possessing abundant elements for light, stress, and hormonal responses. In contrast, Group VI's members show a specialized profile with a strong stress-responsive capacity but relatively fewer light or hormonal elements. Such differential cis-regulatory configurations suggest that although some MsABF genes function as general integrators of environmental cues, others may have evolved more specialized roles in specific stress adaptation pathways.

Figure 5. 

The cis-acting elements of ABF gene promoters in M. sativa. (a) Phylogenetic tree of the MsABF genes. (b) Distribution of the positions of cis-acting elements. (c) Quantitative heat map of cis-acting elements. (d) Bar chart of classified cis-acting elements.

Expression profiling of MsABF genes under various conditions

Transcriptomic analysis revealed distinct spatiotemporal and stress-responsive expression patterns among the MsABF family members. All 37 identified genes showed detectable expression levels across the six examined tissues (root, elongating stem, pre-elongated stem, leaf, flower, and nodule), indicating their broad functional roles in alfalfa's development. However, specific members displayed pronounced tissue-specific expression: MsABF1 and MsABF4 showed significantly higher expression in elongating stems, whereas MsABF18, MsABF20, and MsABF25 were predominantly expressed in the leaves, suggesting specialized roles in stem elongation and leaf development, respectively for the two groups.

Under abiotic stress conditions, the MsABF family exhibited coordinated yet differential regulatory responses (Fig. 6). Following salt stress, eight MsABF genes exhibited marked upregulation. Notably, the same eight genes were also significantly induced under drought stress, indicating a shared molecular response to osmotic stress. Under cold stress, a broader set of 12 genes showed significant induction. Importantly, the 8 salt/drought-responsive genes formed a core subset within the 12 cold-responsive genes, suggesting that these members function as central regulators in general stress adaptation, whereas additional cold-specific regulators (MsABF1, MsABF4, MsABF36, and MsABF44) may be involved in temperature-specific signaling pathways. These expression patterns reveal a hierarchical regulatory architecture: A core set of MsABF genes responds to multiple abiotic stresses, whereas other members exhibit more specialized induction patterns, reflecting both conserved and differentiated functions within this transcription factor family.

Figure 6. 

(a) Expression of MsABF genes in different tissues. (b) Expression of ABF genes in M. sativa under salt stress. (c) Expression of ABF genes in M. sativa under drought stress. (d) Expression of ABF genes in M. sativa under cold stress. CK represents the untreated control. For salt and drought treatments, S1, S3, S6, S12, and S24 and M1, M3, M6, M12, and M24 indicate samples collected at 1, 3, 6, 12, and 24 hours after treatment, respectively. For cold treatment, C2, C6, C24, and C48 indicate samples collected at 2, 6, 24, and 48 hours after treatment.

Protein–protein interaction network analysis

To explore the functional associations of MsABF proteins, we constructed the predicted PPI network using the STRING database and visualized it in Cytoscape (Fig. 7). The resulting network consisted of 35 nodes (including MsABF members and their interaction partners) connected by 99 edges, forming a complex interconnected architecture. Analysis of the network's topology revealed distinct functional modules. Several MsABF proteins, notably MsABF17 and MsABF25 (highlighted in orange), exhibited extensive interactions with M. truncatula orthologs (indicated by cyan or nodes with the prefix 'MEDTR'), suggesting conserved functional complexes within legumes. A distinct subnetwork was formed by MsABF24, which specifically interacts with G7KE11-MEDTR, G7LIV7-MEDTR, and G7KSG0-MEDTR, indicating a potential specialized functional module. According to the sequence homology with their respective M. truncatula orthologs, these hub proteins were speculated to function as molecular scaffolds that mediate the interaction between MsABF transcription factors and downstream stress-responsive signaling components, thereby facilitating the integration and transmission of abiotic stress signals. This suggested that these hub proteins might coordinate crosstalk between distinct MsABF-mediated regulatory pathways.

Figure 7. 

Protein–protein interaction network map. The orange marker is MsABF protein, and the purple node is recognized as a highly intermediate central protein.

Centrality analysis identified key topological hubs within the network. Proteins in the central region (highlighted in purple), such as AOA072VGL4 and AOA072UYK5, ranked highest in betweenness centrality, indicating their critical role as connectors bridging different functional modules. This suggests that these hub proteins may coordinate crosstalk between distinct MsABF-mediated regulatory pathways. The PPI network reveals both conserved interaction patterns (through MEDTR orthologs) and specialized MsABF-specific modules, reflecting the evolutionary conservation and functional diversification of this transcription factor family in regulating stress responses.

Validation of stress-responsive MsABF gene expression by RT-qPCR

To validate the stress-responsive expression patterns identified through transcriptomic analysis, we selected three candidate MsABF genes (MsABF18, MsABF20, and MsABF25) that showed significant responses to all three abiotic stresses and moderate basal expression across tissues for RT-qPCR analysis under controlled stress treatments.

The expression profiles of the three genes exhibited consistent induction patterns across different stress conditions, though with distinct temporal dynamics (Fig. 8). Under drought stress (15% PEG), all three genes showed peak induction at 6 hours post-treatment (hpt), with expression levels approximately 17-, 18-, and 28-fold higher than the control for MsABF18, MsABF20, and MsABF25, respectively. Expression subsequently declined at 12 and 24 hpt. Under cold stress (4 °C), maximum induction was observed at 24 hpt, with six-, five-, and ninefold increases for MsABF18, MsABF20, and MsABF25, respectively. Notably, expression at 6 hpt was higher than at 12 hpt, suggesting a biphasic response pattern.

Figure 8. 

RT-qPCR analysis of MsABF's expression under drought, cold, and salt stress conditions. Significance levels: * p < 0.05; ** p < 0.01, *** p < 0.001, **** p < 0.0001. (a) The expression of MsABF18 under three different stresses. (b) The expression of MsABF20 under three different stresses. (c) The expression of MsABF25 under three different stresses.

In response to salt stress (200 mM NaCl), peak expression occurred at 6 hpt (four-, three-, and fivefold increases), followed by a decrease at 12 hpt and a secondary upregulation at 24 hpt. These RT-qPCR results confirm that MsABF18, MsABF20, and MsABF25 are genuine early-responsive genes under osmotic stress (drought and salt), but they show more sustained induction under cold stress. The differential temporal expression patterns suggest distinct regulatory mechanisms for these genes in response to different abiotic stresses.

In this study, we identified 46 MsABF genes in the alfalfa genome, which were phylogenetically classified into seven distinct groups. The expansion of this family appears to be driven primarily by segmental duplication events, as shown by the 53 duplication pairs identified. This expansion pattern aligns with observations in other plant species, where whole-genome and segmental duplications have been major forces shaping transcription factor families' evolution^[18,19]. Notably, MsABF2, MsABF18, and MsABF20 participated in six segmental duplication events each, suggesting that these genes may have undergone positive selection during alfalfa's evolution, possibly because of their functional importance in stress adaptation. Among them, MsABF18 and MsABF20 exhibited significant differential expression under drought, salt, and cold stress conditions, as determined by both transcriptomic profiling and RT-qPCR validation. Their high duplication frequency, together with their conserved stress-responsive expression patterns, implies that these two genes have been preferentially retained and functionally optimized during evolution to serve as core regulatory nodes in abiotic stress signaling networks^[20]. In contrast, MsABF2 showed no significant stress induction, suggesting it may have undergone subfunctionalization or acquired divergent regulatory roles following duplication. However, previous functional studies have demonstrated that MsABF2 is actively involved in stress responses in alfalfa. It has been reported that MsABF2 directly activates MsLEA-D34 by binding to the ABRE element in its promoter, thereby regulating abiotic stress responses^[21]. Additionally, MsABF2 directly activates MsMYB741 through ABRE element binding, contributing to aluminum stress tolerance be regulating flavonoid biosynthesis^[22]. Furthermore, MsABF2 functions as an upstream regulator of MsDHN1 in an ABA-dependent manner, enhancing aluminum tolerance by promoting oxalate exudation from the root tips^[21]. These findings indicate that although MsABF2 may not be transcriptionally induced under the stress conditions examined in this study, it likely plays important roles in stress adaptation through post-transcriptional regulation or in response to specific stresses such as aluminum toxicity.

All MsABF proteins contain the canonical bZIP domain, confirming their classification within this transcription factor family. The identification of 10 conserved motifs with varying distribution patterns across different members suggests functional diversification. Particularly noteworthy is the extreme variation in motif composition between MsABF29 (containing only one motif) and MsABF9 (containing all 10 motifs). This pattern is consistent with findings in other plant species, where gene family expansion often leads to subfunctionalization or neofunctionalization through structural variations^[23,24]. The conservation of gene structure, with all members lacking UTRs and maintaining relatively stable intron–exon organization, further supports the notion that functional diversification occurred primarily through sequence evolution rather than major structural rearrangements.

Promoter analysis revealed a diverse landscape of regulatory elements within the MsABF family. In total, 826 cis-elements were identified and categorized into four major functional groups with distinct phylogenetic distributions. The concentration of light-responsive elements in Groups I and III, stress-responsive elements in Groups I and VI, and hormone-responsive elements in Groups I and III suggests specialized regulatory strategies among different subfamilies. This differential cis-regulatory architecture likely enables precise spatiotemporal control of MsABF's expression in response to various environmental cues. Similar regulatory complexity has been observed in ABF homologs from other species, where the promoters' architecture reflects their functional specialization^[4]. The predicted PPI network reveals potential functional relationships among MsABF proteins. The extensive interactions between MsABF17, MsABF24, and M. truncatula orthologs suggest conserved protein complexes within the legumes^[25].

Our transcriptomic analysis revealed distinct spatiotemporal expression patterns, with all 46 genes detected across six examined tissues. The tissue-specific expression of MsABF1 and MsABF4 in elongating stems and MsABF18, MsABF20, and MsABF25 in the leaves suggests specialized developmental functions consistent with findings in other plant species, where ABF members regulate organ-specific developmental processes^[11,26].

The RT-qPCR validation of selected stress-responsive genes provides important confirmation of their functional relevance. MsABF18, MsABF20, and MsABF25 showed consistent induction patterns under three different abiotic stresses, though with distinct temporal dynamics. Under drought stress, all three genes peaked at 6 hpt (17-, 18-, and 28-fold increases, respectively), followed by a gradual decline. This rapid response pattern aligns with the known role of ABF/AREB transcription factors as early responders to osmotic stress^[27,28]. Under cold stress, maximum induction occurred at 24 hpt (six-, five-, and ninefold increases respectively), with a biphasic response pattern observed at earlier time points. This delayed response suggests potential involvement in cold acclimation processes, consistent with reports of ABF members participating in temperature stress responses^[29,30]. Salt stress induced different temporal patterns, with peak expression at 6 hpt (four-, three-, and fivefold increases), followed by a decrease at 12 hpt and secondary upregulation at 24 hpt. This biphasic response may reflect distinct regulatory mechanisms for immediate ion homeostasis and longer-term adaptation, similar to observations in other species^[31,32].

This study presents the first comprehensive genome-wide characterization of the ABF/AREB transcription factor family in alfalfa, identifying 46 MsABF genes classified into seven phylogenetic groups. Gene duplication analysis revealed segmental duplications as the primary mechanism of expansion, with 53 duplication pairs identified. Structural analyses showed conserved bZIP domains alongside variable motif compositions, suggesting functional diversification. Promoter analysis identified 826 cis-regulatory elements exhibiting distinct group-specific distributions.

Transcriptomic profiling revealed broad tissue expression with organ-specific patterns, whereas RT-qPCR validation confirmed that selected members (MsABF18, MsABF20, and MsABF25) respond significantly to multiple abiotic stresses, exhibiting distinct temporal dynamics. The identification of these stress-responsive genes indicates important candidates for understanding the mechanisms of stress adaptation in alfalfa.