Abstract

Background

   Alfalfa (Medicago sativa L.) is an important high-quality forage crop.
   Low temperature is an abiotic stress factor that affects the
   distribution and productivity of alfalfa. To further understand the
   molecular response to low temperature, and to identify additional genes
   and metabolic pathways associated with cold tolerance in alfalfa, in
   this study we conducted transcriptome sequencing, weighted gene
   co-expression network analysis, KEGG pathway enrichment analysis, and
   quantitative real-time PCR validation in alfalfa cultivars subjected to
   low-temperature treatment.

Results

   Weighted gene co-expression network analysis revealed that three gene
   modules were significantly negatively correlated with the semi-lethal
   temperature for alfalfa. Genes in the three modules were used to
   construct gene co-expression networks, from which MS.gene46105,
   MS.gene044087, MS.gene76894, MS.gene44620, MS.gene22005, MS.gene045060,
   MS.gene31405, and MS.gene74761 were selected as important genes
   associated with cold tolerance. Quantitative real-time PCR analysis of
   these eight genes validated the reliability of the transcriptome
   sequencing data. In addition, further analysis of the genes within the
   three modules revealed that several transcription factors (AP2/ERF,
   bZIP, C3H, NAC, and others) and metabolic pathways (N-glycan
   biosynthesis, citrate cycle, glycolysis/gluconeogenesis, and carbon
   metabolism, and others) responded well to the low temperature.

Conclusions

   Three gene modules, eight genes, several transcription factors and
   multiple metabolic pathways associated with cold tolerance were
   screened. This results will provide a valuable reference for further
   clarification of the cold tolerance mechanism and breeding for cold
   tolerance in alfalfa.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s12870-024-05987-5.

   Keywords: Alfalfa, Low temperature, Cold tolerance, Molecular response,
   Weighted gene co-expression network analysis

Background

   Low temperature is an abiotic stress factor that limits the
   distribution and productivity of crops globally. The harm caused by low
   temperature to plants includes chilling injury and freezing damage
   [[38]1]. In recent years, abnormal meteorological conditions, such as
   abrupt cooling and extreme cold temperature, have increased in
   frequency. Given the uncertainty associated with their occurrence, low
   temperatures pose a serious threat and challenge to agricultural
   production. Alfalfa, as a high-quality perennial leguminous forage
   crop, plays an important role in livestock production worldwide owing
   to its adaptability, high yield, and nutritional value [[39]2].
   However, alfalfa cultivation is often affected by low-temperature
   stress. Understanding the impact of low temperature on alfalfa and its
   response to cold is important for the formulation of agricultural
   management measures and the breeding of cold-tolerant cultivars.

   Exposure to non-freezing low temperatures for an extended period, a
   process termed cold acclimation, enhances the freezing tolerance of
   alfalfa [[40]3]. During this process, extensive changes in gene
   expression, biochemistry, physiology, and metabolism are observed in
   alfalfa. Many genes associated with cold tolerance have been
   identified, including transcription factors such as DREB/CBF
   (Dehydration-responsive element-binding /C-repeat binding factors),
   CAMTA (Calmodulin-binding transcription activator), and MYB-like family
   members [[41]4–[42]5], as well as genes such as CML (Calmodulin-like)
   [[43]6], ATGs (Autophagy–related genes) [[44]7], eEF-2 (Eukaryotic
   translation elongation factor-2) [[45]8], RLKs (Receptor-like kinases)
   [[46]9], MIPS1 (Myo-Inositol phosphate synthase 1) [[47]10], GolS1
   (Galactinol synthase) [[48]11], INT-like (Inositol transporter-like
   protein) [[49]12], and TIL1 (Temperature-induced lipocalins 1)
   [[50]13]. The proteins or enzymes encoded by genes enhance the cold
   tolerance of alfalfa by accumulation of various metabolites. Osmotic
   substances and cryoprotectants, such as soluble sugars (sucrose,
   raffinose, and stachyose) [[51]14–[52]15] and nitrogen-containing
   compounds (proline) [[53]16], accumulate during cold acclimation. The
   antioxidant system is activated to scavenge reactive oxygen species
   (ROS) that accumulate as a result of cold stress and thus alleviate
   ROS-induced cellular damage [[54]17].

   RNA-sequencing has been used extensively in recent years to study the
   molecular mechanisms of developmental processes and abiotic stress
   responses in alfalfa. A reasonably large number of genes has been
   identified that play critical roles in the response to cold stress.
   Many previous studies, however, were restricted to studying only one or
   two alfalfa cultivars; as a result, the relevant conclusions cannot be
   generalized to lend an overall understanding of the molecular
   mechanisms of cold tolerance in alfalfa [[55]18–[56]20]. In recent
   years, with the accumulation of gene expression data, network-based
   approaches have been employed to explore gene functions or molecular
   mechanisms at the biological process level. Notably, weighted gene
   co-expression network analysis (WGCNA), currently the most efficient
   method, has been widely used to identify correlations among genes in
   large-scale studies. Furthermore, WGCNA assists in characterizing the
   molecular functions of gene modules comprising co-expressed genes. The
   approach involves conducting correlation analysis between co-expressed
   genes and biological traits, effectively identifying hub genes
   associated with the biological traits [[57]21]. WGCNA has been proven
   to be effective for detection of co-expressed modules and hub genes,
   microRNAs, and long non-coding RNAs [[58]22].

   In this study, three alfalfa cultivars were subjected to
   low-temperature treatment to compare their tolerance with semi-lethal
   temperature (LT[50]), followed by transcriptome sequencing of crown
   samples collected at different time points. WGCNA was used to analyze
   the correlation between the LT[50] and transcript abundances, and to
   screen gene modules and hub genes within modules that were strongly
   correlated with cold tolerance. These genes were subjected to
   bioinformatic analysis and quantitative real-time PCR (qRT-PCR)
   validation. In addition, for the genes within the strongly correlated
   gene modules, Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment
   analysis was performed to identify the pathways significantly
   associated with cold tolerance. The findings provide a foundation for
   further clarification of the mechanism of cold tolerance in alfalfa.

Methods

Plant material

   Given the strong relationship between cold tolerance and fall dormancy
   in alfalfa [[59]23], three Chinese alfalfa cultivars that differ in
   fall dormancy rating were selected as experimental materials for
   low-temperature treatment, namely, ‘Xinjiang Daye’ (XJ; semi-dormant;
   identified by Jichun Min et al. in 1990), ‘Zhaodong’ (ZD; dormant;
   identified by Diankui Wang et al. in 1989), and ‘Caoyuan No. 2’ (CY;
   extremely dormant; identified by Yongfu Wu et al. in 1987). These three
   alfalfas are cultivated widely in China, and the seeds we used were
   bought from Jiuquan Daye Seed Industry Co., Ltd. (Jiuquan city, Gansu
   province, China).

Experimental design

   In an experimental field located at Taigu, Shanxi, China, one-year-old
   seedlings of uniform size were selected for each cultivar and
   transplanted to polyvinyl chloride pipes (diameter 10 cm, height
   10 cm). After 1 month, the pipes were moved to a growth chamber
   (LRG-LED-350-D, Hangzhou Lvbo, China) for adaptive growth for 1 week at
   room temperature. The conditions for adaptive growth were set at
   20 °C/16°C (light/dark) with light intensity of 20,000 lx and a
   12 h/12 h (light/dark) photoperiod. For low-temperature treatment, the
   temperature regime was changed to 8 °C/4°C (light/dark), while the
   light conditions were unchanged. During low-temperature treatment,
   crown samples were collected at six time points (0, 12, 24, 48, 96, and
   192 h), with three biological replicates for each cultivar at each time
   point.

Semi-lethal temperature

   The LT[50] was used to characterize the cold tolerance of the alfalfa
   cultivars following the method described by Xu et al. [[60]24]. First,
   the samples were subjected to a gradient freezing process. The freezing
   test was conducted using a constant-temperature circulator (ZX-5 C,
   Shanghai Zhixin, China) at nine gradually decreasing temperatures.
   Subsequently, the relative osmotic permeability of the cellular fluid
   at each frozen temperature was determined by calculating the
   electrolyte conductivity. A logistic equation was used to simulate the
   gradient freezing temperatures and corresponding relative osmotic
   permeability of the cellular fluid. Finally, the LT[50] value was
   determined as the temperature at which the relative osmotic
   permeability of the cellular fluid attained 50%. For the freezing test,
   the nine gradually decreasing temperatures were as follows in
   Table [61]1.

Table 1.

   Nine decreasing temperatures for samples freezing test
   Order of samples freezing test
   Samples for freezing test 1th 2nd 3rd 4th 5th 6th 7th 8th 9th
   0 h 0 -2 -4 -5 -6 -8 -10 -12 -14
   12 h -2 -5 -7 -9 -10 -12 -14 -16 -19
   24 h -2 -5 -7 -9 -10 -12 -14 -16 -19
   48 h -4 -7 -10 -12 -14 -16 -18 -21 -24
   96 h -4 -6 -8 -10 -11 -12 -14 -16 -18
   192 h -2 -4 -6 -8 -10 -12 -14 -16 -18
   [62]Open in a new tab

   Note: the unit of temperatures in freezing test is ℃

Statistical analysis of the semi-lethal temperature

   Statistical analysis of differences in semi-lethal temperatures and
   gene expression levels was conducted using IBM SPSS Statistics for
   Windows, version 19 software (IBM Corp., Armonk, NY, USA). Significance
   was assessed in multiple comparisons using Fisher’s least significant
   difference test, with the significance threshold set to P < 0.05.

Transcriptome sequencing

   Total RNA was extracted from the crown samples of 3 cultivars, 6 time
   points and 3 biological replicates, then mRNA libraries were
   constructed for each sample and sequenced using an Nova Seq 6000
   platform (Illumina, USA). After filtering the raw data, checking the
   sequencing error rate, and estimating the GC content, the clean reads
   were aligned to the alfalfa reference genome [[63]25] using HISAT2
   v4.8.2. The gene alignment was quantified using Feature Counts v1.6.2.
   The unmapped reads were spliced using StringTie [[64]26], and a new
   transcript was assembled and named the novel gene. The new transcript
   is compared with the functional database and plant transcription factor
   database, including Kyoto Encyclopedia of Genes and Genomes (KEGG)
   database ([65]http://www.genome.jp/kegg/), Gene Ontology (GO) database
   ([66]http://geneontology.org/), NCBI nonredundant protein (NR) database
   ([67]http://www.ncbi.nlm.nih.gov), Swiss-Prot (A manually annotated and
   reviewed protein sequence database, [68]http://www.uniprot.org/),
   Translation from European Molecular Biology Laboratory (TrEMBL)
   database ([69]http://www.uniprot.org/), Clusters of Orthologous Groups
   of proteins (KOG) database ([70]http://www.genome.jp/kegg/ko.html),
   Plant Transcription Factor Database (PlnTFDB) database
   ([71]http://plntfdb.bio.uni-potsdam.de/v3.0/) [[72]27] and Plant
   Transcription Factor Database (PlantTFDB) database
   ([73]https://planttfdb.gao-lab.org/) [[74]28], to predict their
   function. The transcript expression levels were expressed as fragments
   per kilobase of transcript per million fragments mapped (FPKM).

Weighted gene co-expression network analysis

   The FPKM values derived from the transcriptome sequence data were
   initially filtered using the Var Filter function of the gene filter
   package for R [[75]29] to eliminate genes with low expression levels
   and genes exhibiting stable expression patterns across all samples. The
   filtered FPKM data and LT[50] values were subjected to correlation
   analysis using WGCNA v1.69 software [[76]30] to identify gene modules
   that were significantly associated with LT[50]. In this analysis, the
   mRNAs and LT[50] values comprised the data for all three alfalfa
   cultivars (CY, ZD, and XJ) at the six sampling time points, without
   distinguishing the results based on the cultivar. Hub genes (crucial
   cold tolerance-related genes) were selected based on the connectivity
   among genes within the modules, and the steps of hub genes selection in
   the three modules are shown in Supplementary Tables [77]S1, [78]S2,
   [79]S3. In order to better understand the effect of hub genes on cold
   tolerance of alfalfa, sequence alignment for the hub genes was
   performed against the NCBI non-redundant protein sequences (NR)
   database to predict the identity of the encoded proteins.

KEGG pathway enrichment analysis

   After screening the gene modules, a KEGG pathway analysis was performed
   on all genes within the gene modules significantly associated with cold
   tolerance to identify the relevant enriched metabolic pathways.

qRT-PCR validation of gene expression

   WGCNA was conducted between the transcriptome data and the LT[50] of
   each sample to screen for genes associated with cold tolerance. After
   identifying cold tolerance-related genes, qRT-PCR was employed to
   validate the expression of the selected genes and the RNA-Seq data,
   three additional alfalfa cultivars, namely, ‘WL168HQ’ (fall
   dormancy = 2), ‘WL343HQ’ (fall dormancy = 4), and ‘WL440HQ’ (fall
   dormancy = 6), were subjected to the same low-temperature treatment as
   described in experimental design. The seeds of three alfalfa were
   provided by Jiuquan Daye Seed Industry Co., Ltd. (Jiuquan city, Gansu
   province, China). Crowns were sampled at 0, 24, and 168 h time points,
   and five biological replicates of each of the samples were prepared for
   qRT-PCR analysis. MiniBEST Plant RNA Extraction Kit (TaKaRa, Dalian,
   China) was used to isolate the total RNA of the whole crown samples.
   Subsequently, a Spectra Max iD^3 multi-mode microplate reader
   (Molecular Devices, CA, USA) was applied to determine the RNA
   concentration of each sample. The first-strand cDNA of each sample was
   synthesized using a PrimeScript™ RT reagent Kit with gDNA Eraser
   (TaKaRa, Dalian, China) according to the manufacturer’s instructions.
   The expression of the cold tolerance-related genes were explored by
   using a CFX Connect™ Real-Time System (Bio-Rad, Los Angeles, CA, USA)
   with TB Green^® Premix Ex Taq™ II (TaKaRa, Dalian, China). Primer
   Premier 6 software (Premier Biosoft International, Palo Alto, CA, USA)
   was applied to design the specific primers for qRT-PCR, and the primer
   sequences are listed in Supplementary Table [80]S4. Ms-Actin2 gene was
   used as the reference gene. The relative quantification method
   (2^−ΔΔCT) was used to calculate the relative expression levels of
   selected genes [[81]31].

Results

Semi-lethal temperature of alfalfa crowns

   Compared with the 0 h time point, the LT[50] of XJ and ZD showed no
   significant change at 96 h, whereas the LT[50] of CY decreased
   significantly (P < 0.05) (Fig. [82]1). Therefore, the ability of CY to
   tolerate cold changed more rapidly than that of ZD and XJ under
   exposure to low temperature. In present study, the minimum LT[50] of ZD
   (− 12.69 °C) was significantly smaller than that of XJ (− 8.85 °C), but
   the minimum LT[50] of CY (− 8.55 °C) was not smaller than that of ZD
   and XJ, and CY might require longer treatment time to reach its minimum
   LT[50]. The LT[50] of XJ and ZD decreased significantly at 192 h
   (P < 0.05), with the LT[50] of XJ reaching − 8.85 °C and that of ZD
   reaching − 12.69 °C. Thus, the cold tolerance of ZD changed more
   rapidly than that of XJ. In addition, the LT[50] of CY reached its
   minimum value at 96 h (− 8.55 °C), and thereafter slightly increased at
   192 h but not significantly (P > 0.05).

Fig. 1.

   [83]Fig. 1
   [84]Open in a new tab

   Semi-lethal temperature (LT[50]) of alfalfa crowns. A smaller LT[50]
   value indicates stronger cold tolerance. Different letters indicate a
   significant difference in LT[50] between different sampling time points

Sequencing and sequence alignment of alfalfa transcriptome

   RNA-sequencing was performed for three alfalfa cultivars subjected to
   cold-stress treatment. The constructed RNA libraries were sequenced for
   54 samples, and a total of 2,398,779,208 raw reads and 2,315,478,432
   clean reads were generated (Supplementary Table [85]S5). Approximately
   90% of the clean reads were mapped to the reference genome
   (Supplementary Table [86]S6). In total, 166,982 unigenes, including
   44,756 novel genes, were assembled (Supplementary Table [87]S7).

Screening of gene modules associated with cold tolerance

   A total of 33,397 genes (Supplementary Table [88]S8) were filtered and
   used for WGCNA. The WGCNA software was employed to calculate the
   correlation coefficients between any two genes, with a correlation
   coefficient threshold of 0.85 considered to indicate a strong
   association. In addition, the weighted values of the gene–gene
   correlation coefficients were analyzed by determining the β power
   value, which strengthens strong correlations and weakens weak
   correlations, leading to a network in which the connections between
   genes followed a scale-free network distribution. The most suitable β
   value selected for this study was 7 (Fig. [89]2A). A hierarchical
   clustering dendrogram was constructed based on the pairwise gene
   correlation coefficients and gene modules were delineated
   (Fig. [90]2B), resulting in the division of the genes for analysis into
   28 modules (Fig. [91]2C). Among these modules, the brown (r = − 0.53,
   P = 3.8 × 10^− 5), darkgreen (r = − 0.6, P = 1.6 × 10^− 6), and red
   (r = − 0.45, P = 6.4 × 10^− 4) modules showed a highly significant
   negative correlation with the cold tolerance of alfalfa (P < 0.01).
   These three modules were selected as the main focus of further
   analyses.

Fig. 2.

   [92]Fig. 2
   [93]Open in a new tab

   Results of weighted gene co-expression network analysis. (A) Selection
   of the soft threshold. The left panel depicts the log k value on the
   vertical axis, representing the square of the correlation coefficient,
   and the horizontal axis is the weight parameter β (soft threshold). The
   right panel shows the mean connectivity of all genes in the
   corresponding gene module on the vertical axis. (B) Hierarchical
   clustering dendrogram, where each color in the tree represents a
   module, and each gene belonging to the same module is indicated by the
   corresponding color. The vertical distance is the distance between two
   genes, and the horizontal distance is not meaningful. (C) Correlations
   between gene modules and the cold tolerance of alfalfa. Darker colors
   indicate stronger correlations. The numbers outside the parentheses are
   the correlation coefficients, and the numbers inside the parentheses
   represent the significance of the correlations

Transcription factors and metabolic pathways response to low temperatures

   Genes included in the brown (2956 genes), darkgreen (186 genes), and
   red (1765 genes) cold tolerance-associated modules (Supplementary Table
   [94]S9) were subjected to a screening of transcription factors and a
   KEGG enrichment analysis. Among all the genes in the three modules, 195
   genes related to transcription factors were identified. Among them, the
   main transcription factors involved are AP2/ERF, bZIP, C3H, NAC, MYB,
   C2H2, and so on. Figure [95]3 shows the types of transcription factors
   involved.

Fig. 3.

   [96]Fig. 3
   [97]Open in a new tab

   Classification of related transcription factors in three gene modules.
   The numbers outside parentheses are the number of related genes
   involved. The numbers in parentheses are the proportion of the number
   of genes (195) in this statistic

   The pathways that involved the genes were classified in five
   categories: cellular processes, environmental information processing,
   genetic information processing, metabolism, and organismal systems. The
   top 20 significantly enriched metabolic pathways in the three gene
   modules (Fig. [98]4A–C) included N-glycan biosynthesis, citrate cycle,
   glycolysis/gluconeogenesis, carbon metabolism, oxidative
   phosphorylation, proteasome, peroxisome, ribosome, arginine and proline
   metabolism, and MAPK signaling pathway – plant.

Fig. 4.

   [99]Fig. 4
   [100]Open in a new tab

   Scatterplots of KEGG pathway enrichment analysis of genes within three
   cold tolerance-associated gene modules. (A) Brown, (B) darkgreen, and
   (C) red gene modules. The vertical axis represents KEGG pathways, and
   the horizontal axis is the ratio of ‘the number of differentially
   expressed genes annotated with the pathway term’ to ‘the total number
   of genes annotated with the pathway term’ (rich factor). A larger rich
   factor indicates a higher degree of enrichment. The size of the points
   corresponds to the number of genes enriched in the pathway. The color
   of the points represents the significance of pathway enrichment,
   assessed using the q-value. A q-value closer to 0 is indicated by a
   redder color, representing a more significant enrichment. The figure
   displays the top 20 significantly enriched pathways

Screening of crucial genes associated with cold tolerance

   A gene co-expression network was constructed based on the strength of
   connections among the genes within the three modules most significantly
   associated with cold tolerance, namely, the brown, darkgreen, and red
   modules (Fig. [101]5A–C). The eight hub genes with highest connection
   strength were selected as important genes associated with cold
   tolerance in alfalfa, comprising one gene from the brown module, three
   genes from the darkgreen module, and four genes from the red module.
   The sequences of the eight selected genes were queried against the NCBI
   NR database to predict the identity of the encoded protein
   (Supplementary Table [102]S10).

Fig. 5.

   [103]Fig. 5
   [104]Open in a new tab

   Gene co-expression networks of partial genes within the (A) red, (B)
   darkgreen, and (C) brown gene modules. The red nodes represent the hub
   genes, and the yellow nodes represent other genes in the gene networks
   with relatively high connection strength

Expression of eight important genes associated with cold tolerance

   Among the eight genes associated with cold tolerance, the expression
   levels of all genes except MS.gene76894 exhibited a similar trend
   throughout the experimental period, with the expression level
   increasing with prolonged duration of low-temperature treatment
   (Fig. [105]6A–H). The expression levels of MS.gene46105, MS.gene44620,
   MS.gene22005, MS.gene045060, MS.gene31405, and MS.gene74761 increased
   rapidly, and a significant increase in expression level at 192 h
   compared with that at 0 h was observed (P < 0.05). In contrast,
   MS.gene044087 exhibited non-significant changes in expression levels
   during the 0–96 h treatment period, with a slow rate of change
   observed, but a significant increase (P < 0.05) in expression level was
   detected in the 96–192 h treatment period. The gene MS.gene76894 did
   not show a significant change in expression level throughout the
   experimental period (P > 0.05).

Fig. 6.

   [106]Fig. 6
   [107]Open in a new tab

   FPKM values derived from transcriptome sequencing data for eight cold
   tolerance-associated genes at six time points during cold treatment.
   The x-axis is the sampling time points and the y-axis is the gene
   expression level. Different letters indicate a significant difference
   between time points within the same cultivar (P < 0.05)

qRT-PCR validation of eight important genes associated with cold tolerance

   A qRT-PCR analysis was conducted to validate the expression patterns of
   the eight selected cold tolerance-associated genes in three different
   cultivars of alfalfa, namely WL168HQ, WL343HQ, and WL440HQ
   (Fig. [108]7A–H). Except for MS.gene76984, the remaining seven genes
   exhibited a gradual increase in expression level in all three or
   certain of the cultivars during the treatment period. This trend was
   essentially consistent with the change in FPKM values derived from the
   transcriptome sequencing data (Fig. [109]6A–H). The gene MS.gene76984
   showed a gradual increase in expression level in WL168HQ, whereas an
   initial increase followed by a decline was observed in WL343HQ and
   WL440HQ (P < 0.05). The genes MS.gene46105 and MS.gene44620 showed a
   significant increase in expression level in WL168HQ and WL343HQ
   (P < 0.05), but no significant change in expression level was observed
   in WL440HQ. The gene MS.gene22005 exhibited a significant increase in
   expression level in WL168HQ (P < 0.05), but no significant change in
   expression level was observed in WL343HQ and WL440HQ. Thus, among the
   responses to cold treatment, more pronounced changes in the expression
   of cold tolerance-related genes were observed in alfalfa cultivars with
   weaker dormancy and stronger cold tolerance, such as WL168HQ.

Fig. 7.

   [110]Fig. 7
   [111]Open in a new tab

   Expression levels of eight selected cold tolerance-associated genes
   during cold treatment in the alfalfa cultivars ‘WL168HQ’, ‘WL343HQ’,
   and ‘WL440HQ’. The x-axis is the sampling time points and the y-axis is
   the gene expression level. Different letters indicate a significant
   difference in gene expression level between different time points
   within the same cultivar (P < 0.05)

Discussion

Effects of identified genes and predicted proteins on cold tolerance of
alfalfa

   In present study, the cold tolerance-associated gene MS.gene46105 was
   predicted to be a member of the glycoprotein family. In addition, KEGG
   enrichment analysis in our study revealed significant enrichment of
   metabolic pathways, such as N-glycan biosynthesis, various types of
   N-glycan biosynthesis, types of O-glycan biosynthesis, and protein
   processing in the endoplasmic reticulum (ER). In plant cells,
   glycosylation is among the most studied post-translational events. It
   can be of two types, N- or O-glycosylation [[112]32]. Of these two
   types, the N linkage is the most common [[113]33]. Typically, most
   proteins located in the extracellular space and membrane systems of
   plant cells undergo N-glycosylation [[114]34]. For proteins,
   N-glycosylation has an important impact on glycoprotein structure,
   conformation, stability, solubility, immunogenicity, and enzyme
   activity. The N-glycosylation of proteins plays functional roles in
   regulating endoplasmic reticulum quality control (ERQC), stress
   tolerance, signal recognition, and signal transduction [[115]35].
   Glycoproteins are both assembled and modified within the endomembrane
   system, i.e., the ER and the Golgi apparatus, before transport to their
   ultimate location within or outside the cell [[116]36]. ERQC, through
   the binding of N-glycoproteins to the ER chaperones, recognizes
   unfolded or misfolded protein intermediates to be retained in the ER
   for an additional cycle of refolding. If refolding fails, the
   accumulated misfolded proteins are subsequently targeted and
   translocated for ER-associated degradation. The correctly folded
   glycoproteins are transported to the Golgi apparatus, plasma membrane,
   and other subcellular compartments, where complex-, hybrid- and
   paucimannose-type N-glycans are formed in mature proteins [[117]33,
   [118]37]. In the research on Arabidopsis, N-glycan degradation of
   stress-regulated proteins in the ER is observed responding to cold
   treatment. Thus, proper regulation of ERQC to trigger N-glycan
   degradation via the unfolded protein response is required to maintain
   the correct folding of stress-responsive proteins, which may be an
   important mechanism for plant adaptation to cold stress [[119]35].
   Under exposure to cold stress, in the basal portion of the leaf sheath
   in rice, changes are observed in the glycosylation and phosphorylation
   profiles of calreticulin, a crucial protein that regulates the quality
   control of other proteins [[120]34]. Thus, ER stress responses are
   activated by misfolded proteins that accumulate in the ER under abiotic
   stress, which is one crucial mechanism for plants to mitigate
   stress-induced damage and confer stress tolerance [[121]38].

   Regarding O-glycoproteins, arabinogalactan proteins and extensions are
   two O-glycosylated superfamilies of hydroxyproline-rich glycoproteins
   that have received considerable attention in recent years owing to
   their diverse functions associated with plant development, and
   responses to biotic and abiotic stresses [[122]36]. In contrast to
   N-glycoproteins, the functions of O-glycoproteins in plant responses to
   abiotic stress are poorly understood. Based on the above analysis, we
   hypothesize that, under cold stress, alfalfa likely undergoes
   glycosylation of various stress-responsive proteins in the ER to ensure
   their correct folding, maintain their proper three-dimensional
   structure, protect their functionality, enhance the stability of the
   endomembrane system, and ultimately improve cellular cold tolerance.

   In NR database, MS.gene044087 was predicted to be SNAK1. SNAK1, a
   gibberellin-regulated family protein, was initially isolated from
   potato tubers as a cysteine-rich antimicrobial peptide, which is an
   ancient and crucial component of innate immunity [[123]39–[124]40].
   Given the homology of its amino acid sequence, SNAK1 is classified as a
   member of the Snakin/Gibberellic Acid Stimulated in Arabidopsis (GASA)
   protein family [[125]39]. Silencing of SNAK1 leads to a significant
   reduction in plant height, decrease in leaf size, and severe changes in
   leaf morphology in potato [[126]41]. In addition, the ROS content
   increases in the leaves, the ascorbic acid (AsA) content significantly
   decreases, and salicylic acid (SA) content increases [[127]42]. AsA is
   a common non-enzymatic antioxidant in plant cells that has the ability
   to scavenge ROS and participates in plant responses to abiotic stress
   [[128]43–[129]44]. SA, a phenolic compound, plays a role in signal
   transduction under both biotic and abiotic stress, and can induce the
   production of plant defense metabolites, including antioxidant
   components [[130]45]. Furthermore, SA can enhance plant cold tolerance
   by activating antioxidant enzymes, such as ascorbate peroxidase and
   glutathione reductase, to eliminate ROS generated by cold stress
   [[131]46]. Thus, the silencing of SNAK1 indirectly reduces the cold
   tolerance of potato. Whether MS.gene044087 can enhance cold tolerance
   by activating the antioxidant defense system, further validation is
   necessary to determine.

   The gene MS.gene76894 is highly likely to be a LON peptidase N-terminal
   domain and RING (Really Interesting New Gene) [[132]47] finger protein
   3. LON peptidases are ubiquitous, multidomain, ATP-dependent enzymes
   with both highly specific and nonspecific protein binding, unfolding,
   and degrading activities [[133]48]. The ATP-dependent LON peptidases
   are evolutionarily conserved in all living organisms and catalyze rapid
   turnover of short-lived regulatory proteins and many damaged or
   denatured proteins [[134]49]. In Escherichia coli, LON peptidases are
   responsible for degrading most newly synthesized abnormal proteins
   [[135]50]. The N-terminal domain of the LON peptidase is essential for
   substrate (target protein) selection and recognition during its
   functional activity [[136]51–[137]53], and the N-terminal region of LON
   family members is not highly conserved [[138]50]. The RING finger is a
   type of zinc finger protein. The RING finger domain can be categorized
   into classical and modified RING finger subfamilies based on the
   combination of cysteine and histidine residues, with the C3HC4 type
   belonging to the classical RING finger subfamily [[139]54]. Regarding
   the function of the RING domain, the promoter of AtATL15, a C3HC4-type
   RING zinc finger gene, in Arabidopsis is responsive to AsA [[140]55].
   In Brassica napus, BrRZFP1, a C3HC4-type RING zinc finger protein,
   strongly induces abscisic acid (ABA), and overexpression of BrRZFP1 in
   tobacco increases plant tolerance to cold, salt, and dehydration stress
   [[141]56]. ABA plays a central role in seed dormancy, abscission, and
   transmission of abiotic stress signals [[142]57]. ABA regulates many
   critical processes in plant growth, development, and adaptation to
   biotic and abiotic stress, such as stomatal regulation and the
   expression of defense-related genes [[143]58]. In the present study,
   MS.gene76894 was strongly associated with cold tolerance of alfalfa.
   However, further validation is required to determine whether the gene
   affects cold tolerance by regulating the synthesis of phytohormones,
   such as ABA, through the RING finger domain.

   The gene MS.gene44620 was predicted to be a glycine-rich RNA-binding
   protein 3 (GR-RBP). Different types of RBPs play roles in the
   regulation of RNA metabolism in various cellular processes [[144]59].
   Conserved RNA-binding motifs include RNA recognition motif,
   glycine-rich motif, arginine-glycine-glycine, RGG box motif, zinc
   finger motif, and double-stranded RNA-binding motif [[145]60]. The
   typical feature of GR-RBPs is the presence of a glycine-rich domain at
   the C-terminus [[146]61]. A glycine-rich zinc finger RNA-binding
   protein in rice, OsRZ2, exhibits RNA chaperone activity under cold
   stress [[147]59]. The glycine-rich RNA-binding protein NbGRP7 regulates
   the function of the Gpa2/Rx1 protein at the post-transcriptional level
   and acts as a pro-immune component in effector-triggered immunity
   [[148]62]. Increasing evidence suggests that up-regulation of
   glycine-rich RNA-binding proteins participates in the mediation of
   stress responses under various stresses, such as low temperature, ABA,
   SA, injury, ultraviolet radiation, salt, pathogen infection, and water
   deficit [[149]61, [150]63]. The gene MS.gene44620 may function as a
   member of the GR-RBP family and play a crucial role in enhancing the
   cold tolerance of alfalfa by mediating stress responses under cold
   stress.

   Tryptophan synthase is composed of two subunits, α and β [[151]64]. The
   α subunit cleaves indole-3-glycerol phosphate into indole and
   glyceraldehyde-3-phosphate, whereas the β subunit catalyzes the
   condensation of indole and serine to produce tryptophan [[152]65]. In
   this study, MS.gene22005 was predicted to encode tryptophan synthase
   beta chain 1. Tryptophan synthase beta subunit 1 plays a crucial
   coordinating role in plant growth and stress responses by maintaining
   the homeostasis of tryptophan and ABA [[153]66]. Melatonin, a
   derivative of tryptophan, controls hydrogen peroxide accumulation in
   plants by scavenging excess ROS, and by enhancing antioxidant enzyme
   activities and ascorbate–glutathione recycling [[154]67]. In addition,
   melatonin enhances plant tolerance to various stresses, such as cold,
   heat, salt, and drought, by regulating genes in the DREB/CBF, HSF, SOS,
   and ABA pathways [[155]67–[156]68]. Expression of MS.gene22005 may be
   involved in tryptophan synthesis in alfalfa or contribute to melatonin
   production, thereby influencing the cold tolerance of alfalfa through
   participating in a series of physiological and biochemical responses.

   The cold tolerance-associated gene MS.gene045060 aligned with protein
   disulfide isomerase (PDI) in the NCBI database. PDI is primarily
   localized in the ER of eukaryotic cells and is capable of catalyzing
   the formation and isomerization of disulfide bonds in substrate
   proteins [[157]69–[158]70]. It prevents the aggregation of unfolded or
   partially folded protein precursors, and serves as both a
   thiol-disulfide oxidoreductase and a molecular chaperone in many
   organisms [[159]71–[160]72]. Typical PDIs contain four thioredoxin-like
   domains, a, a′, b, and b′. The enzyme activity of PDI is directly
   associated with the CXXC motif in the a and a′ domains, and mutations
   in the conserved cysteines of this motif cause the loss of enzymatic
   activity [[161]69, [162]73]. PDI proteins in higher plants are involved
   in signaling pathways and play a role in responses to abiotic stresses,
   such as salt, drought, and low temperature [[163]69].

   Based on the present gene co-expression networks, MS.gene31405 and
   MS.gene74761 were also identified as potential cold
   tolerance-associated genes in alfalfa. Although there is currently
   limited research on the two genes, they warrant attention in future
   studies of the cold tolerance mechanism in alfalfa.

Effects of related transcription factor on cold tolerance of alfalfa

   The three modules, brown, red and darkgreen, involve 14 classes of
   transcription factors. These transcription factors could regulate plant
   adaptation to stress by regulating gene expression, substance
   metabolism or hormone signaling pathways. For instance, AP2/ERF
   (APETALA2/ETHYLENE RESPONSE FACTOR) is involved in transcriptional
   regulation under abiotic stress, as well as in both hormone-dependent
   and hormone-independent signaling pathways [[164]74]. AP2/ERF also
   participates in responses to abiotic stress mediated by auxin,
   gibberellin, and cytokinin [[165]75]. The bZIP (Leucine zipper)
   transcription factor family plays a multifaceted role in plants’
   response to abiotic stress. For example, the abscisic acid response
   element binding factor is a member of the bZIP transcription factor
   family and plays a central regulatory role in ABA signaling pathway
   [[166]76]. Furthermore, bZIP transcription factors play crucial roles
   in multiple biological processes such as, gibberellic acid signaling,
   metal ion signaling, pathogen perception, and sucrose signaling
   [[167]77]. Zinc-finger proteins, a superfamily of proteins with a
   typical structural domain that coordinates a zinc ion and binds nucleic
   acids, participate in plant development and tolerance to salt, drought,
   flood low temperature and oxidative stresses caused by non-biological
   factors [[168]78]. Most zinc fingers are C2H2-type or CCCC-type, and
   C3H is a type with a CCCH structure [[169]79]. In three modules of
   present research, 10 and 12 genes are involved in C2H2 and C3H,
   respectively. N-acetylcysteine (NAC) transcription factors can bind to
   target gene DNA-binding domains to activate or inhibit the expression
   of marker gene. This has been shown to be an important tool for
   enhancing plant responses to water stress/drought stress across
   different plants [[170]80].

Effects of metabolites and metabolic pathways on cold tolerance of alfalfa

   We identified three gene modules strongly associated with cold
   tolerance in alfalfa. KEGG enrichment analysis of all genes in the
   three modules revealed numerous metabolic pathways potentially
   associated with cold tolerance, such as N-glycan biosynthesis and
   carbon metabolism (Fig. [171]4A–C). These metabolic pathways included
   several involved in energy synthesis and metabolism, such as the
   citrate cycle, glycolysis/gluconeogenesis, carbon metabolism, oxidative
   phosphorylation, and pyruvate metabolism. Previous research on cold
   acclimation in alfalfa has shown that, compared with cultivars with
   lower cold tolerance, hardy cultivars exhibit a more pronounced
   increase in root functions (respiration, nodulation capacity, and
   various enzyme activities) [[172]81–[173]82]. Consistent with the
   present findings, genes involved in the ATP synthase complex and
   ribosome biogenesis are likely to be involved in the alfalfa response
   to cold stress [[174]19]. Thus, it is highly probable that energy
   metabolism is strengthened to cope with low-temperature stress. The
   capacity of energy metabolism under low temperatures may serve as a
   basis for evaluation of the cold tolerance of alfalfa accessions.

   Under stress, plants produce excessive ROS, prompting the initiation of
   antioxidant defense systems to eliminate ROS and protect cells from
   oxidative damage. Notably, the peroxisome metabolism pathway was
   significantly enriched in the present study. Peroxisomes contain a
   variety of enzymes, primarily oxidases, catalases, and peroxidases,
   which, through their combined activity, can mitigate the harmful
   effects of toxic substances, such as hydrogen peroxide, on cells
   [[175]83]. Thus, under low-temperature stress, peroxisomes may reduce
   cell damage caused by hydrogen peroxide, thereby protecting cells and
   enhancing cold tolerance in alfalfa.

   The mitogen-activated protein kinase (MAPK) cascade consists of
   different combinations of at least three protein kinases, MAPKKK
   (MAP3K/MEKK/MKKK), MAPKK (MKK/MEK), and MAPK (MPK), which activate each
   other in a sequential manner through phosphorylation [[176]84]. Cold
   stress induces the expression of C-repeat-binding factors (CBFs) for
   which Inducer of CBF expression 1 (ICE1) is the main regulatory factor
   for CBFs. Some MPKs can affect the stability of ICE1, thereby
   influencing plant cold tolerance [[177]85]. Specific members of the
   MAPK family in plants are reported to be activated by ABA and play
   roles in protection of cell signaling, disease resistance and abiotic
   stress tolerance, and regulation of catalase 1 in ROS homeostasis
   [[178]84]. Thus, the MAPK cascade signaling pathway likely plays a
   crucial role in protection against low-temperature stress in alfalfa.

   Genes from the three cold tolerance-associated modules were
   significantly enriched in many other important metabolic pathways, such
   as arginine and proline metabolism, biosynthesis of amino acids
   metabolism, arachidonic acid, fatty acid metabolism, fatty acid
   biosynthesis, carbon fixation in photosynthetic organisms, and fructose
   and mannose metabolism. Many previous studies have investigated the
   effects of amino acids, fatty acids, and carbohydrates on the cold
   tolerance of alfalfa [[179]3, [180]86–[181]90], to which the reader is
   referred for further details. Further research is warranted on their
   involvement in cold tolerance mechanisms in alfalfa.

Conclusion

   In this study, three gene modules strongly associated with cold
   tolerance of alfalfa were identified. In three modules, eight crucial
   genes associated with cold tolerance were identified, and they were
   MS.gene46105, MS.gene044087, MS.gene76894, MS.gene44620, MS.gene22005,
   MS.gene045060, MS.gene31405, and MS.gene74761. Upon conducting
   additional statistical analyses on the three gene modules, we
   identified 14 classes of transcription factors, encompassing AP2/ERF
   (Apetala2/Ethylene response factor), bZIP (Leucine zipper), C3H
   (Zinc-finger protein), and NAC (N-acetylcysteine) families. These
   factors are pivotal in the plant’s response to abiotic stress,
   including resistance to low temperatures. Subsequent KEGG enrichment
   analysis elucidated the involvement of multiple metabolic pathways,
   including N-glycan biosynthesis, citrate cycle,
   glycolysis/gluconeogenesis, and carbon metabolism, in the cold response
   of alfalfa.

Electronic supplementary material

   Below is the link to the electronic supplementary material.
   [182]Supplementary Material 1^ (5.1MB, xlsx)
   [183]Supplementary Material 2^ (325KB, xlsx)
   [184]Supplementary Material 3^ (3.6MB, xlsx)
   [185]Supplementary Material 4^ (9.3KB, xlsx)
   [186]Supplementary Material 5^ (10.2KB, xlsx)
   [187]Supplementary Material 6^ (11.5KB, xlsx)
   [188]Supplementary Material 7^ (59.8MB, xlsx)
   [189]Supplementary Material 8^ (16.4MB, xlsx)
   [190]Supplementary Material 9^ (2.7MB, xlsx)
   [191]Supplementary Material 10^ (9.3KB, xlsx)

Author contributions

   HX conceived and designed the study. ZZ performed most of the
   experiments. HX and ZZ wrote and revised the manuscript. QZ, YG, YX,
   and JC assisted with data acquisition. YL and XH assisted with data
   analysis.

Funding

   This study was financially supported by the National Natural Science
   Foundation of China (No. 32101450), the Award Fund for Excellent
   Researchers Work in Shanxi Province (No. SXBYKY2021017), and the
   Science and Technology Innovation Foundation of Shanxi Agricultural
   University (No. 2020BQ64). The funding body had no contribution to the
   design of the study, to the collection, analysis, and interpretation of
   data or to the writing of the manuscript.

Data availability

   The RNA-seq raw data that support the findings of this study have been
   deposited in the NCBI SRA database with the primary accession code
   PRJNA1090639.

Declarations

Ethics approval and consent to participate

   Plant materials used in this study were cultivated on the experimental
   plot of Shanxi Agricultural University by our research team, and all
   the members have no objection to this use.

Consent for publication

   Not applicable.

Competing interests

   The authors declare no competing interests.

Footnotes

   Publisher’s note

   Springer Nature remains neutral with regard to jurisdictional claims in
   published maps and institutional affiliations.

   Hongyu Xu and Zipei Zhang contributed equally to this work.

References