Differential expression analysis of genes in the root crown of alfalfa under low-temperature stress


 Background: Alfalfa ( Medicago sativa ) is a perennial forage crop widely cultivated in northern China. The root crown is an important storage organ of alfalfa, especially in the wintering process, as it is closely related to winter hardiness. However, the molecular mechanism underlying the winter hardiness of the alfalfa root crown remains unclear. To investigate these gaps in knowledge, the RNA sequencing (RNA-Seq) technology was used to identify critical genes related to winter hardiness. Results: In this study, the winter survival rate of the Lomgmu 806 variety was approximately 3.68-fold higher than that of the Sardi variety. We sequenced the transcriptomes of the root crown of the two alfalfa varieties. Among the 57,712 unigenes identified, 2,299 differentially expressed genes (DEGs) were upregulated, and 2,143 DEGs were downregulated. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotations showed that 1,159 unigenes were mainly annotated in 116 pathways. Seven DEGs belonging to the “plant hormone signalling transduction” pathway, the “peroxisome” pathway and transcription factor family (MYB, B3, AP2/ERF, and WRKY) and involved in alfalfa winter hardiness were identified. As a result, the expression patterns of seven DEGs were verified by real-time quantitative PCR (RT-qPCR) analyses, which verified the reliability of the RNA-Seq analyses. Conclusions: The RNA-Seq data revealed the gene regulation response of alfalfa to low-temperature stress, which provides a valuable resource for the further identification and functional analysis of candidate genes related to winter hardiness in alfalfa. Furthermore, these data provide references for future in-depth studies of winter hardiness mechanisms in alfalfa.

4 winter hardiness mechanisms in alfalfa.

Background
Alfalfa (Medicago sativa L.) is one of the most popular perennial legume forages and is known as the "king of the forages", due to its high yield, rich nutritional value and palatability. This forage is not only the most widely cultivated in the world [1,2] but also the most widely distributed in China [3,4]. Alfalfa varieties introduced from abroad have the advantages of good quality and high yield, but their winter survival rates or spring returning green rates are low, which seriously limits their industrial development in northern China [2,4]. In fact, there are also problems in the cultivation of alfalfa varieties that cannot safely survive under harsh winter conditions in the midwestern United States, Canada, Italy, Russia and other countries [5]. Due to the lack of winter-hardy varieties in some cold regions, understanding the molecular mechanism of low-temperature stress signal transduction pathways is very important for breeders aiming to improve the production of alfalfa.
Plants can survive safely in winter and enhance freezing resistance through a lowtemperature adaptation process [6]. Previous studies on the mechanisms of plant responses to low temperature have mainly focused on plant physiology, phytohormones, and changes in low temperature-induced gene expression by regulating transcription factors (TFs) [7]. Plants undergo various changes under lowtemperature stress, and the key response to low-temperature stress depends on the activation of molecular networks involved in stress signal transduction pathways and the expression of specific sets of stress-related genes [8]. At present, the most widely studied C-repeat-binding factor (CBF) signalling pathway including CBF1, 5 CBF2, and CBF3, and CBF TFs regulates the expression of COLD-RESPONSIVE (COR) functional genes, which are components of the low-temperature signal transduction pathway, thereby improving the cold tolerance of plants [9,10]. A number of studies have shown that the CBF signalling pathway can improve the lowtemperature tolerance of various plants, such as Arabidopsis thaliana [11,12], Triticum aestivum [13], Glycine max [14], and Malus domestica [15]. The overexpression of MtCBF3 cloned from Medicago truncatula induces COR gene expression and enhances the low-temperature tolerance of transgenic plants [16].
Recent studies have shown that CBFs may play an important role in regulating the low-temperature tolerance of M. sativa [17,18]. This suggests that the regulation of the CBF pathway may have a potential role in improving the low-temperature tolerance of alfalfa, although it may not be the only method by which lowtemperature tolerance is improved.
The plant hormone signal transduction pathway plays a central role in the plant lowtemperature stress response and controls the CBF-dependent and CBF-independent pathways [19]. For example, abscisic acid (ABA) is an important abiotic stressregulating hormone that plays an important role in abiotic stress signalling.
Furthermore, the previous assumption suggested that ABA did not affect the expression of the CBF gene. Therefore, it has been suggested that the lowtemperature response via ABA control is independent of CBF regulation [20].
However, a recent study suggested that ABA may affect the expression of the COR genes by regulating CBF transcription and may play an important role in plant lowtemperature stress in plants [21]. Moreover, some functional genes also play an important role in protecting plants from cold and/or freezing stress. These genes are involved in signal transduction and transcriptional regulation. For example, TFs play 6 a functional role by regulating the expression or status of stress-related genes. At present, some TFs have been identified, such as members of the MYB, B3, AP2/ERF, and WRKY families, and serve as important regulators of abiotic stress responses in plants [22][23][24][25]. In addition, in the plant response to abiotic stress, reactive oxygen species (ROS), which are harmful to cell membranes, proteins, and biological macromolecules, accumulate in plant cells destroying the homeostasis of the organism and causing serious damage to plants [26]. Therefore, the exact roles of the response factors of alfalfa in the wintering period, including transcriptional regulation, the plant hormone signal transduction pathway, and ROS scavenging, remain to be elucidated, and their functions should be studied in detail.
Recently, transcriptome sequencing has become essential for functional gene annotation, novel gene discovery, differential gene expression, and molecular marker research [27][28][29]. RNA sequencing (RNA-Seq) enables the study of gene expression at the transcriptome level and identifies genes involved in plant-specific biological processes [30][31][32]. For example, with the introduction of a new generation of RNA-Seq, transcriptome sequencing technology has been used to identify lowtemperature stress response genes in Ipomoea batatas [33], Populus tomentosa [34], Capsicum annuum [35], and Magnolia wufengensis [36]. Moreover, some recent studies have reported that the use of transcriptome sequencing technology to analyze the expression of genes in the seedlings [7], taproots [37], and crown buds [38] of alfalfa, which provides valuable resources for functional genomics research on alfalfa cold tolerance. However, the understanding of genetic regulation information regarding low-temperature stress remains limited. Therefore, in this study, alfalfa root crowns were collected during the wintering period for RNA-Seq to obtain genetic information on the response of alfalfa to low-temperature stress.
Moreover, the specific molecular mechanisms underlying the low-temperature resistance of the alfalfa root crown during the wintering period remain unclear.
Therefore, the systematic analysis of the low-temperature resistance of the Longmu 806 and Sardi varieties during winter not only aided in the understanding of the response mechanism of the alfalfa root crown during wintering but also provided a reference for functional genomics research on low-temperature stress in alfalfa.

Winter survival
In the spring of 2019, Longmu 806 alfalfa exhibited the larger returning green plants than Sardi alfalfa in the experimental field (Fig. 1a). In September 2018, there was a significant difference in plant height (PH) between Longmu 806 and Sardi. The PH of Sardi alfalfa was significantly higher than that of Longmu 806 alfalfa (P<0.05; Fig. 1b). The winter survival rate of Longmu 806 (98.09%) was approximately 3.68-fold higher than that of Sardi (P<0.05, Fig. 1c).

De novo transcriptome assembly
Total RNA extracted from the root crown of the Longmu 806 and Sardi varieties during the wintering stage were constructed into six libraries for high-throughput sequencing. To ensure the quality of the information analysis, the original data were filtered. Among the raw reads, low-quality reads, reads containing adapter sequences, and reads containing low-quality bases were discarded, and 459,865,666 clean reads with a total of 68,515,729,607 nucleotides (nt) were generated from the six sequencing libraries. A total of 114,567 unigenes with an N50 value of 1,092 nt were assembled. The maximum length, mean length and minimum length were 15,687 nt, 740 nt and 201 nt, respectively (Fig. 2).

Identification and analysis of differentially expressed genes (DEGs)
With the Nr annotation results, BLAST2GO (2.5.0) software was used to perform GO functional annotations of unigenes. A total of 4,442 DEGs were assigned to one or more ontologies by the standard of |log 2 fold change| >1 and P-value < 0.05 (Additional file 1: Table S1; Fig. 4). There were 2,299 upregulated unigenes, and the other 2,143 unigenes were downregulated (Longmu 806 vs Sardi). These DEGs were used for the subsequent analyses.

9
To better understand the function of genes differentially expressed in alfalfa during the wintering period, GO functions were used to classify the function of DEGs (Additional file 2: Table S2). A total of 4,442 DEGs were summarized in three main functional categories: biological process, cellular component, and molecular function (Fig. 5). In the biological process category, "cellular process" (GO:0009987), "metabolic process" (GO:0008152) and "biological regulation"

Stress response of TFs to low-temperature stress
A total of 34 TF families containing 1,364 unigenes were identified (Additional file 4,

Validation of DEGs data by real-time quantitative PCR (RT-qPCR) analysis
To verify the reliability of the transcriptome sequencing data, differences in the expression of ten genes were detected using RT-qPCR, and these genes were well characterized by the NCBI Nr database (Fig. 7a)

Discussion
Low-temperature stress is one of the main abiotic stresses affecting plant growth and development and alfalfa yield [2]. Generally, low-temperature freezing damage is likely to occur under severe winter conditions, which could result in reduced forage yields in the following year and seriously affect forage production. Therefore, the ability for regrowth in early spring following severe winter conditions reflects the degree of cold resistance in forage [40]. In this study, we also found that Longmu 806 alfalfa had more returning green plants than Sardi alfalfa in the early spring (Fig. 1a). In addition, a recent study has indicated that the winter survival rate of a fall dormant alfalfa (Maverick) was markedly higher than that of a nondormant alfalfa (CUF101), suggesting that fall dormant alfalfa has a stronger cold tolerance than non-dormant alfalfa [37]. Similar results were obtained in this study.
The winter survival rate of Longmu 806 (fall dormant alfalfa) was significantly higher than that of Sardi (non-dormant alfalfa; Fig. 1c). In summary, the winter hardiness of Longmu 806 was stronger than that of Sardi. Because the winter hardiness of the alfalfa variety is a complex trait, it is not only affected by the growth environment but also related to the comprehensive genetic regulation of the variety [41]. To our knowledge, previous studies have reported the use of RNA-Seq technology to analyze DEGs related to alfalfa leaves, taproots, and crown buds [37,38,42]. However, unlike previous studies, to enhance our understanding of the comprehensive genetic regulation mechanism of alfalfa low-temperature resistance in alfalfa, we focused on the expression patterns of low-temperature responsive genes in the root crown of alfalfa.
Hormones in alfalfa in response to low-temperature stress Hormones are the initiation factors of winter hardiness gene expression, and hormone-mediated abiotic stress responses involve multiple mechanisms. In perennial plants, the phytohormones ABA, ethylene, JA and auxin play important roles in regulating plant growth to adapt to low-temperature stress [43]. ABA is an important stress-regulating hormone in plants under low-temperature stress. A recent study found that the levels of ABA increased in plants treated with low temperature, and ABA-dependent and ABA-independent pathways can regulate lowtemperature-responsive genes [44]. SnRK2 is a positive regulator of the ABA 13 signalling pathway, and the role of the ABA-activated SnRK2 protein kinase in lowtemperature stress signalling has been reported in Kandelia obovata [45]. In this study, the RNA-Seq data showed that the SnRK2 gene was significantly upregulated in the Longmu 806 variety but downregulated in the Sardi variety. This result suggested that Longmu 806 was more likely to accumulate ABA than Sardi. The ABAbinding receptor inactivates PP2C and PP2C-mediated dephosphorylation of the SnRK2 kinase in Arabidopsis, which phosphorylates the basic leucine zipper (bZIP) TFs (ABFs/AREBs).. ABFs bind to ABA-responsive promoter elements (ABREs) to activate the transcription of downstream ABA-responsive genes [46,47]. Therefore, the overexpression of the SnRK2 gene in Longmu 806 regulated the low-temperature stress response by inducing the expression of ABA-responsive genes. Similar results have also been found with the overexpression of TaSnRK2.3 in Arabidopsis, which results in an improved root-system structure and significantly enhanced the tolerance of this species to freezing stress [48]. In addition, low-temperature stress and ABA synergistically (in Vitis vinifera buds) increased the expression of CBF/DREB1 TFs (VvCBF2, VvCBF3, VvCBF4, and VvCBF6) [49], and increased the cold resistance of the Longmu 806 variety, indicating that Longmu 806 is more likely to resist complex low-temperature environments than Sardi.
EIN3 is involved in ethylene signal transduction, and ethylene signalling promotes the transcription of multiple ethylene response factor (ERF) genes, ultimately guiding growth and physiological responses to abiotic stress [50]. The RNA-Seq data revealed that the EIN3 gene was differentially expressed between the Longmu 806 and Sardi varieties. The RT-qPCR analysis showed that the expression level of EIN3 in the Longmu 806 variety was approximately 2.2-fold higher than that of Sardi variety. This indicated that Longmu 806 was more likely to accumulate ethylene 14 than Sardi. EIN3 promotes the transcription of a variety of ERF genes that ultimately regulate the expression of downstream stress-responsive genes. For example, overexpressing the CBF/DREB genes in tobacco and tomato enhanced freezing tolerance by regulating the expression of the COR genes [51]. In addition, a recent study revealed several mutants (etr1-1, ein4-1, ein2-5, ein3-1, and ein3 eil1) in the ethylene signalling pathway, in which the ein3 eil1 mutant induced the expression of CBF genes (CBF1, CBF2, and CBF3) under low-temperature stress and then enhanced the freezing tolerance of Arabidopsis [52]. Interestingly, high concentrations of ABA inhibit root growth in Arabidopsis by increasing ethylene accumulation [53]. Therefore, although the mechanisms underlying the complex cross talk between ABA, ethylene, and low-temperature signalling remain unclear, it appears likely that the upregulation of ABA-responsive genes also helps the low- tolerance [54]. The RNA-Seq data showed that the JAZ gene in Longmu 806 alfalfa was upregulated but the downregulation of the gene in Sardi alfalfa was not significant. This indicated that Longmu 806 alfalfa may reduce JAZ expression through ICE-CBF regulatory pathways, releasing CBF TFs to induce the expression of the COR genes and ultimately improving alfalfa cold resistance [54]. Similarly, the expression of the ARF gene in Longmu 806 was upregulated compared with that in 15 Sardi. Shu et al. [55] also found that the transcription of the ARF gene in M.sativa was upregulated under cold and freezing stresses, which indicated that the conserved miRNAs might regulate the alfalfa low-temperature stress response by controlling the alfalfa developmental process. Although the ARF gene is not directly involved in the alfalfa root overwintering, it may affect the elongation and dormancy of alfalfa roots under low-temperature stress by regulating auxin [38].
Therefore, although the ARF and JAZ genes of Longmu 806 were upregulated, their exact roles under low-temperature stress remain to be elucidated, and their functions require further study.
The response of DEGs to low-temperature stress revealed the relationship between the "plant hormone signal transduction" pathway (ko04075), "MAPK signalling pathway-plant" pathway (ko04016), and winter hardiness. The DEGs associated with ABA and ethylene biosynthetic pathways (including SnRK2 and EIN3) in Longmu 806 alfalfa were upregulated. Hormones, as signalling molecules, play very important roles in regulating CBF gene expression in response to cold and freezing stresses in M. sativa (cv. Zhaodong) [18]. In summary, it is likely that Longmu 806 alfalfa has stronger winter hardiness than Sardi alfalfa due to the high levels of accumulated plant hormones. Increased expression of the SnRK2 and EIN3 genes caused an increase in the ABA and ethylene contents and then induced the expression of CBF genes, which may be a reason for increased cold resistance in the Longmu 806 variety.
TFs involved in the low-temperature stress response TFs play important roles in the response to abiotic (low temperature, salt, and drought) stresses and directly control the expression of specific sets of stressresponsive genes [56,57]. Most TF families identified play a vital role in plant responses to low-temperature stress, including the MYB, AP2/ERF, B3, and WRKY families [22][23][24][25]. A recent study has reported that many TF subfamilies belong to the AP2/ERF family, such as the AP2 subfamily, DREB subfamily, ERF subfamily, RAV subfamily and others, and a high percentage were involved in the low-temperature (cold and/or freezing) stress response [58]. In this study, we identified that the DEGs related to the GO terms "cellular process", "metabolic process", "cell part", "cell", "organelle", "nucleic acid binding transcription factor activity", and "binding" were enriched. The AP2/EREBP gene of the AP2/ERF family was significantly upregulated in low-temperature stress responses of the Longmu 806 variety, and similar results were obtained in M. truncatula. In particular, the AP2/ERF TF family includes a tandem array of Mt-ERF genes (also called MtCBF genes) on chromosome 6 that functions in the response to cold and freezing stresses in M. truncatula [22]. This result verifies previous reports in Camellia sinensis that five CsAP2/ERF genes from each AP2/ERF family (AP2, ERF, DREB, and RAV subfamilies) responded to temperature stresses [59]. In brief, this result confirmed that AP2/EREBP TFs have a positive response to low-temperature stress.
The AP2/EREBP gene was upregulated in the Longmu 806 variety compared to the Sardi variety, suggesting that Longmu 806 has enhanced cold resistance.
The MYB TFs have been shown to play positive roles in the abiotic stress signalling process [60]. In the MYB family, the JAmyb geneis induced by high salinity, high osmotic stresses and ROS, indicating that JAmyb is responsible for the abiotic stress response [61]. The RT-qPCR analysis showed that the expression level of the TF JAMYB in the Longmu 806 variety was approximately 2.5-fold higher than that in the Sardi variety. This indicated that JAMYB may be a response factor to lowtemperature stress, but its specific function still needs further study. In addition, previous studies have reported that the WRKY family is involved in the lowtemperature stress response. For example, the expression of TFs in most WRKY families (WRKY 30,33,41,46,48,51,53,65,70) is upregulated in Hevea brasiliensis under low-temperature stress (24 h cold treatment at 4 °C) [62]. In particular, the WRKY 65 gene expression in cassava was upregulated under lowtemperature stress [63], which was consistent with the results of our study. The RNA-Seq data showed that the expression of WRKY 65 was upregulated in the Longmu 806 variety but downregulated in the Sardi variety. Increased expression of the WRKY 65 gene can be quickly and significantly induced by ABA and cold treatments [63], indicating that the WRKY 65 gene may positively regulate the cold signalling pathway in an ABA-dependent manner to improve the cold tolerance of Longmu 806 alfalfa. Similar results have also been confirmed in the B3 family. It has been reported that the AP2/ERF and B3 domain-containing TF gene RAV1 was induced by low-temperature stress according to microarray analysis [24,64], and thus may regulate plant growth under low-temperature stress. This indicated that RAV1 may play an important role in restraining root growth under low-temperature stress in the Longmu 806 variety.
Antioxidant defense system-related genes in response to lowtemperature stress Under low-temperature stress, ROS destroy cellular components, disturbing the homeostasis of plants and causing programmed cell death [26]. To protect plants from ROS damage, plants have established a very complex physiological system to scavenge ROS [65]. For example, superoxide dismutase (SOD) is a plant-specific oxidation-reduction enzyme widely present in plants. SOD catalyzes the disproportionation of superoxide anions to oxygen, which protects cells from superoxide poisoning. Increased antioxidant enzyme activity can effectively protect plants against low temperature [66]. In this study, the DEGs associated with SOD1 was upregulated in the "peroxisome" (ko04146) pathway, and the DEGs related to the GO terms "cellular process", "response to stimulus", "catalytic activity", and

Conclusion
In summary, our study identified 4,442 DEGs, and many potential low-temperatureresponsive transcription genes. We found that the antioxidant defense system (SOD1) was involved in the"peroxisome" pathway, and the upregulation of the SOD1 gene may play an important role in improving the low temperature resistance of alfalfa. Moreover, six candidate DEGs were involved in the "plant hormone signal transduction" pathway and in TF families (MYB, B3, AP2/ERF, and WRKY), directly protecting alfalfa from low-temperature stress and increasing the tolerance of alfalfa to a severe winter environment. The cold resistance of alfalfa led to the upregulation of a large number of genes, which may be a protective mechanism to ensure the survival of alfalfa in the winter. This study has improved our understanding of the mechanism underlying alfalfa winter hardiness. In addition, our RNA-Seq data greatly enriched the alfalfa gene resources, which provides a reference for the future in-depth studies of winter hardiness in alfalfa.

Unigene annotation
The unigenes were annotated using BLASTX alignment with an E-value threshold of

Analysis of DEGs
Gene expression analysis of unigenes were calculated and normalized by fragments per kilobase per million reads (FPKM). The DEGs were analyzed using DESeq2 software. After P-value threshold testing using the false discovery rate (FDR), genes with an FDR ≤ 0.05 and a log 2 fold change> 1 were taken as thresholds to select the DEGs. The DEGs were subjected to further GO enrichment analysis and KEGG pathway analysis. TaKaRa, Tokyo, Japan) and an ABI Prism 7500 Sequence Detector (Applied Biosystems). Primers for RT-qPCR were designed using Primer Premier 6.0 software (Premier, USA). All primer sequences are listed in the supplementary material, Table   23 S5. A 20 μl reaction mix was set up containing 50 ng total RNA and 0.4 μM each primer. The RT-qPCR analysis was performed in three technical repeats, and the program was set at 95 °C for 3 minutes, followed by 35   Additional file 2: Table S2. GO functional annotations and the number of differentially expressed genes statistics.

RT-qPCR analysis
Additional file 3: Table S3. KEGG pathway annotation of differential expressed genes in the root crown of two alfalfa varieties.
Additional file 4: Table S4. Transcription factors family statistical table.
Additional file 5: Table S5. Primers used for Real-Time quantitative PCR (RT-qPCR) analysis. Figure 1 The  Expression analysis of the response of ten genes to low-temperature stress by real-time quan 40 Supplementary Files This is a list of supplementary files associated with the primary manuscript. Click to download.