Genome-wide identification and expression analysis of DREB genes in alfalfa (Medicago sativa) in response to cold stress.

Dehydration-responsive element-binding proteins (DREBs) belong to members of the AP2/ERF transcription factor superfamily, which has been reported to involve various abiotic-stress responses and tolerance in plants. However, research on the DREB-family is still limited in alfalfa (Medicago sativa L.), a forage legume cultivated worldwide. The recent genome-sequence release of the alfalfa cultivar "XinJiangDaYe" allowed us to identify 172 DREBs by a multi-step homolog search. The phylogenetic analysis indicated that such MsDREBs could be classified into 5 groups, namely A-1 (56 members), A-2 (39), A-3 (3), A-4 (61) and 13 (A-5 (13), thus adding substantial new members to the DREB-family in alfalfa. Furthermore, a comprehensive survey in silico of conserved motif, gene structure, molecular weight, and isoelectric point (pI) as well as gene expression was conducted. The resulting data showed that, for cold-stress response, 33 differentially expressed MsDREBs were identified with a threshold of Log2-fold > 1, and most of which were transcriptionally upregulated within 48 h during a cold treatment(s). Moreover, the expression profiling of MsDREBs from two ecotypes of alfalfa subspecies i.e. M. sativa ssp. falcata (F56, from a colder region of Central Asia) and M. sativa ssp. sativa (B47, from Near East) revealed that most of the cold-stress responsive MsDREBs exhibited a significantly lower expression in F56, leading to a proposal of the existence of a distinct mechanism(s) for cold tolerance regulated by DREB-related action, which would have been evolved in alfalfa with a genotypic specificity. Additionally, by examining the transcriptome of a freezing-tolerance species (M. sativa cv. Zhaodong), eight DREBs were found to be implicated in a long-term freezing-stress adaptation with a great potential. Taken together, the current genome-wide identification in alfalfa points to the importance of some MsDREBs in the cold-stress response, providing some promising molecular targets to be functionally characterized for the improvement of cold tolerance in crops including alfalfa.

Introduction

Cultivated alfalfa (Medicago sativa L.) represents one of the most important legume forages. Besides their role serving as a livestock feed, such legume plants can also be directly applied as an organic fertilizer to not only provide soil nitrogen (N) derived from symbiotic N-fixation but also to improve soil quality. However, as observed in many other crops, the growth of the cultivated alfalfa suffers often from diverse abiotic stresses, among which a low temperature resulted from e.g. ‘late-spring freezing’ or an abnormally sharp temperature change is a critical factor limiting severely a higher biomass production and quality of alfalfa.[1-3]

To cope with a severe challenge of a detrimental temperature (e.g. more than 5 degree lower than that required for an optimal growth) that may frequently occur in ecosystems, plants have evolved complex and diverse systems to improve their cold-stress resilience, for instance by shifting from a growth to dormant phase or changing morphological characteristics to synchronize with climatic rhythms. At a molecular and/or physiological level, under low temperature or cold conditions, endogenous regulating systems modify the membrane properties e.g. fluidity through alteration of lipid content and composition, misfolding of proteins, increase in anti-oxidative system activities as well as operating other essential strategies to maintain plants in a relatively inactive status.[4,5] In general, such mechanisms are initiated by a cold induced transient increase in a cytosolic concentration of Ca2+, which is mediated by Ca2+-permeable channels and detected by numerous receptors and/or kinase proteins. In addition, the primary cold-stress signal modifies also some second messengers or other signal molecules, such as cyclic nucleotides (cAMP and cGMP), amino acids (e.g. Glu and Met) and various reactive oxygen species (ROS), which may trigger the expression of transcription factors (TFs) required for activation of the cold acclimation regulatory networks at the transcriptome level.[6,7]

The CBF (C-repeat binding factor) and DREB family have been documented to play important roles in response to drought, low-temperature, and high-salt stress in many plant species. Liu et al. reported that DREB1 and DREB2 induced cold- and drought-responsive gene expressions,[8] and this signal transduction pathway is essential for plant response to cold stress.[9] Other studies have indicated that many other TFs (e.g. PIF, MPK, EIN, and CRPK) are also involved in the DREB-dependent signaling processes, emphasizing the importance of DREB proteins in plant cold acclimation.[10] Indeed, transgenic approaches with a manipulation of DREBs’ expression in crops such as tobacco, barley, tomato and Medicago Truncatula can significantly improve plant growth and tolerance to cold stress.[11-14] To date, the putative sequences of DREBs and their AP2/ERF (APETALA 2/ethylene-responsive element) superfamily genes have been identified in many plant species, including Arabidopsis, wild cane, mulberry, wheat and mungbean. For alfalfa plant, although Jin et al.[15] identified 59 DREB/CBF-like proteins, the authors were not able to use more complete genome sequence data, which have been released only recently from the cultivar “Zhongmu No. 1” and “XinJiangDaYe”. As reported for these two cultivars, the 816-Mb high-quality chromosome-level haploid genome sequence (n = 8) and a length of 2738 Mb for 32 (4 n = 32) chromosomes have newly assembled, respectively.[16,17] This may indicate that a large complex genome of autotetraploid alfalfa cultivar would potentially possess more genes due to an evolutional duplication event(s).

In the present study, we identified 172 MsDREB proteins containing the AP2 domain with highly conserved amino acid residues (i.e. 14th valine (V) and 19th glutamic acid (E)) from M. sativa “XinJiangDaYe” genome sequencing dataset. Phylogenetic and conserved motif analysis could classify these genes into 5 subfamilies (A-1 to A-5), but not including A-6 subcluster. Gene structure, conserved motifs, molecular weight, and isoelectric points (pIs) were simultaneously investigated. Furthermore, we analyzed the expression pattern of MsDREB genes in response to cold stress (22°C to 4°C) using RNA-seq datasets from possible websites; 33 core MsDREBs were found to be significantly upregulated during cold treatment, and such MsDREBs were further characterized as low-expression genes in cold tolerance alfalfa cultivar “sativa” but not in “falcata”. Additionally, a transcriptomic analysis with alfalfa cultivar Zhaodong during its long-term cold acclimation unraveled certain core MsDREB genes that may be involved in cold/freezing tolerance.

Materials and methods

Homologs sequence identification and characterization of DREB family in alfalfa

For comprehensive identification and analysis of DREB gene family in alfalfa, sequences of putative DREB genes with 14thV and 19thE were extracted from National Center for Biotechnological Information (NCBI, https://www.ncbi.nlm.nih.gov/) GenBank via BLAST search against Arabidopsis DREBs from the Arabidopsis International Resource (TAIR, https://www.arabidopsis.org/). The draft genome data and annotation information of alfalfa cultivar “XinJiangDaYe” were download from the figshare data repository released on 20–07-2019 (https://figshare.com/projects/whole_genome_sequencing_and_assembly_of_Medicago_sativa/66380). MsDREB genes were first searched by using TBtools “BLAST GUI Wrapper”[17] to obtain candidate DREB-family members. Second, the HMM (Hidden Markov Model) profile PF00847 (AP2 domain) was retrieved from the Pfam database (https://pfam.sanger.ac.uk/) to identify the putative DREB genes from previous BLAST filtered candidates. Basic parameters were predicted using the ProtParam tool (https://web.expasy.org/protparam).[18]

Phylogenetic analysis of DREB proteins

To analyze the sequence feature of DREBs, 208 putative DREB homologs (36 from A. thaliana and 172 from M. sativa) were loaded in the MEGAX version 10.2.6 (https://www.megasoftware.net) using default “ClustalW” parameters for the sequence alignment. Phylogenetic trees were then constructed with default “Neighbor-joining” algorithm parameters, and the best ‘.nwk’ phylogenetic tree was made by iTol (https://itol.embl.de).[19]

Analysis of conserved motifs (MEME), gene structure, and domains (NCBI CDD) with MsDREB genes

Conserved motifs of MsDREBs proteins were predicted using the MEME suite version 5.4.1 (https://meme-suite.org/meme/index.html)[20] with default parameters, except the “motifs should find”, which was set to 10. The coding sequences (CDS) and the structure of all genes were graphically displayed with TBtools function “Gene Structure View”. The gene structure information of “XinJiangDaYe” was downloaded from figshare.com (https://figshare.com/articles/dataset/genome_fasta_sequence_and_annotation_files/12327602?file=22771946).

Expression analysis of MsDREB genes during plant development and in response to cold stress

The raw RNA sequencing data was obtained from NCBI Sequencing Read Archive (SRA, https://www.ncbi.nlm.nih.gov/sra) through the SRA toolkits “prefetch” (version 2.8.0), SRP055547 for development in two alfalfa cultivars, SRP144299 for cold stress and SRP060503.[21-23] Raw data (raw reads) in fastq format were first qualified with FastQC program for Q20, Q30, GC-content and sequence duplication level; the data were then processed in Hisat2 version 2.2.1[24] for read alignment to the genome of alfalfa “XinJiangDaYe”.[17] The reads were subjected to fragments per kilobase of transcript per million fragments mapped (FPKM) conversion to obtain the expression value of genes and transcripts. In-house R scripts were used to analyze gene expression and generate heatmaps. Transcripts (with a fold change (FC) ≥ 2 and an adjusted p-value < 0.05) were identified as differentially expressed genes (DEGs) with DEseq2. A heat map was created using ggplot2, reshape2, and dplyr packages in R version 4.1.2.[25]

Results and discussion

Genome-wide identification and phylogenetic analysis of DREB family genes in M. sativa

To identify the MsDREB genes in the autotetraploid alfalfa M. sativa genome, the BLAST search was performed by using TBtools with the dataset of cultivar “XinJiangDaYe”, query data of known Arabidopsis DREBs with conserved AP2 domain containing 14th valine and 19th glutamic acid residue were processed at an E-value cutoff of 1E-5. Initially, 616 candidate proteins carrying AP2 domain were identified by HMM PF00847; MEGAX alignment was used to verify the conserved AP2 domain and AA residues. A total of 172 MsDREBs with their peptide length ranging from 79 to 505 AA residues were eventually characterized. The calculated molecular weight (MW) of the MsDREBs is in a range of 8.68 kDa to 56.56 kDa, and their isoelectric point (pI) was predicted to be from 4.55 to 10.92 (Table S1). Herein, we have completely surveyed and characterized in M. sativa a large DREB family containing much more members (172) than that of previously reported (59 DREB/CBFs) by Jin et al. (2019).[15]

To explore the phylogenetic relationship of MsDREB genes, 172 MsDREBs and 36 AtDREBs were aligned and processed in MEGAX using ClustalW with proper default parameters (The protein sequences are given in Table S2).[26] Unrooted trees were constructed using the neighbor-joining (NJ) method with specific default settings. The result showed that the MsDREB genes could be separated into 5 groups (except for the A-6 subfamily with distinct sequence motifs, see below), namely A-1 (with 56 members), A-2 (39), A-3 (3), A-4 (61) and A-5 (13) (Figure 1). Concerning again the previously BLAST-predicted 616 MsDREB proteins, we found that most members of A-6 subfamily contained also the AP2 domain, but the conserved 14th V and 19th E were mutated to other amino acids. To date, proteins with the AP2 domain consisting of 14th V and 19th E have been more likely characterized as DREBs, but those with an alanine (A)- and aspartate (D)-containing AP2 motif are assigned to ERF proteins.[27] Because the specific AP2 recognizes and binds to the DRE (dehydration-responsive element) cis-regulatory sequence in regulatory networks, such a AP2 domain represents a typical feature for DREB genes. Thus, in this study, we further used the 172 MsDREB-gene datasets for the preliminary analysis of alfalfa plant in response to cold stress.

Figure 1.

Phylogenetic tree of DREB proteins from Arabidopsis thaliana and Medicago sativa L. Evolutionary relationships were constructed using the neighbor joining (NJ) method (See Materials and Methods). The shaded color indicates the different subfamily.

Conserved motifs and gene structures

To explore the distribution and structural diversification of conserved motifs of MsDREB proteins, we processed a maximum of 10 motifs in silico prediction using MEME online tools and renamed the motifs using a number from 1 to 10 (Figure 2 and Figure S1). Only the A-3 subfamily has all 10 motifs, and the others share parts of the motifs. However, the motif diversity is relatively conserved in all groups, with slight differences among subfamilies. The motif 1 was characterized in all members of MsDREBs. Gene structure analysis showed that the number of exons in MsDREB genes varied from 1 to 4 (Figure 2 and Figure S1). The ‘gene62759’ processes 4 exons and large intron regions. In summary, a large variation of motif numbers and content among different subfamilies of MsDREBs was revealed, different from their gene structures with a less variation, for which 137 out of 172 MsDREB genes had only one exon or without an intron.

Figure 2.

Analysis of the conserved motifs and gene structure of A1-3 subfamily.

Expression patterns of 33 MsDREB genes in response to cold stress

To know if or how the transcription of DREB genes is responsive to cold stress, a publicly available RNA-seq dataset related to time-depend (0 h (CK), 2 h, 6 h, 24 h, 48 h) cold-stress (4°C) treatments of M. sativa was analyzed.[22] The expression pattern of all 172 MsDREBs was obtained from the dataset; only 33 genes were identified as DEGs at a threshold of ‘|Log2fold| > 1’ (p-value < 0.05) using the data from the CK as a reference (Figure 3). Almost all DREBs were continuously induced at the transcriptional level during the tested cold stress, except for gene022547 (showing a downregulation after 24 h and 48 h cold treatment) that was also clustered separately in the phylogenetic tree (Figure 3). The 32 MsDREBs were divided into 3 subclades; genes in the clade I (7 members) were transcriptionally upregulated after 24 h cold treatment, the mRNA level of clade II (17 members) increased initially but relatively declined after 24 h, and the expression of clade III (8 members) was significantly upregulated, with more than 10-fold higher than that of those in the clade I and II (Figure 4). It has been known that DREB family genes are regulated by various types of TFs.[28] In Arabidopsis, CAMTA (calmodulin-binding transcription activator) are Ca2+-dependent trans-activator.[29] Previous studies described that CAMTA3 positively regulated DREB by binding to the cis-element CM2 motif of CBF2 under cold stress and the expression of CBF family genes were significantly downregulated up to approximately 50% in camta1 camta2 double mutant.[30] As a cytosolic oscillation of Ca2+ often reflects an early event in the cold stress response, our results suggest that clade II/III DREB genes would be involved in the rapid action at a low temperature, in which such genes were markedly upregulated during first 2 h cold treatment (Figure 4). Interestingly, the regulation of CBF1/2 expression by CAMTA3/5 is only determined in the case of a sharply decreased temperature (e.g. directly from 22°C to 4°C) but not of gradually lowering of the temperature (e.g. decreasing 3°C every 10 min).[31] Moreover, there is other evidence showing that numerous DREBs play important roles in response to a long-term cold stress, e.g. SsDREB102 and SsDREB105 for E3-mediated ubiquitination during the cold stress in Saccharum spontaneum.[32] In addition, the A-1 subfamily gene CBF2 and DREB displayed differential expression patterns in response to cold; genes in cultivar ‘Chinese Spring’ were only transcriptionally induced in a relative short-term and then followed by rapid decrease in their mRNA abundance;[33] in wheat cultivar “Mironovskaya 808”, the genes were shown their expression at consistently high cold-upregulation levels, similar to that of top half, most clade II DREBs and clade III genes[33,34] (Figure 3). For clade I DREBs, a comparable expression pattern was observed for wheat DREB genes “EVM0034637_A5” and “EVM0026054_A4”; however, their functions as well as associated regulatory networks were largely unclear.[35] Thus, it has been suggested that genetic variations may result from potent differences in the regulatory motifs present in promoters of genes in different cultivars.[36] Although many DREBs have been documented to play diverse and critical roles in plant response to abiotic stimuli including heat, drought and cold,[8,37] for example, the expression of Arabidopsis DREB1/CBF are induced by cold, while the transcription of DREB2 genes are generally upregulated by dehydration, salinity and high temperature,[8,38-40] 33 out of 172 MsDREBs in alfalfa have been characterized here at least as cold-responsive elements, and the rest MsDREBs might participate in other stress regulatory networks, which may be interesting to be explored in the future.

Figure 3.

Clustering approach-based heatmap of MsDREB genes. The heatmap represents the expression pattern of genes in response to cold stress. Relative expression levels were calculated as a Log2-fold change against CK (See Materials and Methods). The red color shows an upregulation of a given gene, and the green indicates a downregulation. Labeled number in each tile is the expression level “FPKM”.

Figure 4.

The expression level of MsDREB genes in different tissues. The expression level is shown by the intensity of blue color. Labeled number in each tile is “FPKM”.

Expression difference between M. sativa subspecies during plant development

To extend our understanding of possible biological functions of the preliminarily selected cold-responsive MsDREB genes, a de novo transcriptome data assembly from two alfalfa subspecies (M. sativa ssp. sativa-B47 from the Middle East, M. sativa ssp. falcata-F56 from colder region Central Asia) was inspected, and RNA-seq data from different tissues were analyzed.[21,41] In this study, we comprehensively examined the expression level of 33 MsDREBs in specific tissues including root, leaf, flower, elongation stem internode, post-elongation stem internode, and nitrogen-fixing nodule (Figure 4). In general, most MsDREB genes showed a low or undetectable expression (FPKM = 1–0) in the cold-tolerant cultivar F56, with a respective number of 31, 25, 25, 17, 21 or 22 MsDREBs in the root, leaf, flower, elongation stem internode, post-elongation stem internode or nodule. Differently, numerous MsDREBs exhibited a relatively high transcription level in some tissues of the cold-sensitive cultivar B47, e.g. gene 85966 and gene 38367 in the leaf and gene 36782 in the root (Figure 4). These results clearly indicate a distinct and ecosystem-dependent expression pattern of 33 MsDREBs in B47 and F56. Being a cold area-adaptive subspecies, since some genes implicate in certain valuable agronomic traits (e.g. winter hardiness, creeping root habit) from F56 plant have already been assigned to cultivated alfalfa,[41] and since such MsDREBs were further characterized in this study as cold-responsive genetic components under a short-term cold stress, we assume that an existence of a cold-induced regulatory network would lead to the different cold response contributed by the action of MsDREBs in different subspecies. However, this hypothesis needs to be experimentally confirmed at a molecular-physiological level. To date, a comparison of gene expression between sensitive and tolerant cultivars has been extensively applied to assess and verify abiotic responsive regulators in many plants, e.g. cowpea, rice, and Chinese jujube under cold stress, wheat and barley subjected to various growth constraints,[42-45] resulting in a successful identification of many potential cold-stress regulatory components. Thus, we expect that those MsDREB genes with differential expression identified between B47 and F56 would be worthy candidates for further functional characterization.

Transcriptome analysis of M. sativa cv. Zhaodong grown adaptively under a long-term freezing condition

M. sativa cv. Zhaodong is a major forage legume cultivar that can survive in the northern China, where temperatures may reach around −30°C in winter. Song et al. (2016) provided a systemic analysis of adaptive mechanisms for frost tolerance using RNAseq assay.[23] M. sativa cv. Zhaodong were cultivated in September and harvested in mid-winter after freezing adaptation during a two-month temperature decrease. In this adaptation process, the antioxidant defense system was believed to confer freezing tolerance in a rapid way, and numerous potential freezing-sensing and signal transduction components were revealed.[23] Here, we focused on the expression pattern of previously identified 33 cold-tolerance related MsDREBs. As shown in Figure 5, the genes in clade 3 were differentially expressed at high transcription levels (FPKM > 10) under the long-term cold/freezing treatment as compared to that of in control samples without cold acclimation. For clade 2, the most genes, whose expressions were previously analyzed to be rapidly upregulated by cold stress but decreased after 24 hours (Figure 3), were not substantially expressed in the freezing acclimation alfalfa (Figure 5). Interestingly, transcripts of the clade 1 genes were highly accumulated in freezing acclimation plants (Figure 5). Such results suggest that the functional expression of clade one-third members would contribute to a long-term freezing fitness of alfalfa plants, while clade 2 genes might be required for an early response (e.g. within 24 h cold) to cold threat. With experimental evidence demonstrating that the overexpression of MtCBF1 (also termed MtDREB1C) triggered by the CaMV 35S RNA promoter leads to a constant and remarkable enhancement of freezing tolerance of Medicago truncatula and Chinese rose,[14] we would like to suggest that the eight MsDREB genes in the clade 3 (Figure 3) would be interesting molecular targets available for the improvement of plant ability to acclimate both cold stress and long-term freezing through genetic manipulation.

Figure 5.

Clustering approach-based heatmap of MsDREB genes in freezing acclimation alfalfa.

In conclusion, here we have genome-widely identified in alfalfa plant 172 putative MsDREB genes using updated whole-genome sequence datasets, detailly characterized their genetic lineage and gene structures as well as coding-proteins that consist of a conserved AP2 domain with 14th valine and 19th glutamic acid residue. In particular, 8 out of 33 differentially expressed MsDREBs in the clade 3 were characterized as potential cold-stress responsive components, which may serve as crucial trans-activators participating in or controlling a signal transduction network(s) regulating growth adaption of crops including alfalfa to cold and/or long-term freezing environments. In the future study, it may be interesting and necessary to functionally appreciate such MsDREB genes at molecular and physiological levels. Thus, we believe that our findings with adding many new members to the DREB-family in alfalfa should provide valuable information for further profoundly understanding of DREBs’ action in plant cold-stress behavior, assisting us to take a proper approach aiming at improving crop cold-tolerance.

Click here for additional data file.