Deciphering the ceRNA Network in Alfalfa: Insights into Cold Stress Tolerance Mechanisms

TL;DR

This study explores how non-coding RNAs and ceRNA networks help alfalfa tolerate cold stress, offering insights into molecular mechanisms and potential targets for improving cold resistance.

Contribution

The study identifies a ceRNA network and key genes involved in cold stress tolerance in alfalfa using RNA-seq and GWAS.

Findings

Ribosome and starch metabolism pathways were significantly enriched with differentially expressed mRNAs under cold stress.

A ceRNA network involving 28 lncRNAs, 8 circRNAs, and 11 miRNAs was constructed to explain cold stress responses in alfalfa.

MS.gene53818 (MsUAM1) was identified as a critical candidate gene for cold stress response through GWAS and RNA-seq analysis.

Abstract

Abiotic stress of cold is one of the limitation factors that hinder the production of alfalfa (Medicago sativa). Although there are a large number of studies suggesting that non-coding RNAs (ncRNAs) play an important role in plant response to abiotic stress, the mechanism by which ncRNAs and competing endogenous RNAs (ceRNAs) influence the low-temperature tolerance of alfalfa remains understudied. In this study, we integrated whole-transcriptome RNA-seq and genome-wide association studies (GWASs) to identify cold stress-related metabolic pathways and candidate genes, differentially expressed (DE) mRNAs, microRNAs (miRNAs), long non-coding RNAs (lncRNAs), and circular RNAs (circRNAs). Degradome sequencing was used to verify the ceRNA network under cold stress. A total of 46,936 DEmRNAs were identified. Ribosome (ko03010), amino sugar and nucleotide sugar metabolism (ko00520), ribosome biogenesis in eukaryotes (ko03008), circadian rhythm–plant (ko00270), and starch and sucrose metabolism (ko00500) were the top five KEGG terms with the highest p-value, enriching the most number of DEmRNAs. MS.gene53818 (MsUAM1) was considered to be the critical candidate gene for alfalfa response to cold stress by conjoint analysis of GWASs and DEmRNAs. A total of 223 DEmiRNAs, 1852 DElncRNAs, and 13 DEcircRNAs were identified under cold stress. Functional analysis indicates that they play important roles in GO terms such as leaf development (GO:0048366), DNA-binding transcription factor activity (GO:0003700), central vacuole (GO:0042807), response to auxin (GO:0009733), and water channel activity (GO:0015250), as well as in KEGG pathways such as plant hormone signal transduction, starch and sucrose metabolism, and flavone and flavonol biosynthesis (ko00944). A ceRNA network comprising 28 DElncRNAs, 8 DEcircRNAs, 11 DEmiRNAs, and 23 DEmRNA triplets was constructed. In this study, mRNAs and ncRNAs were identified that may be involved in alfalfa’s response to cold stress, and a ceRNA regulatory network related to cold stress was established, providing valuable genic resources for further research on the molecular mechanisms underlying alfalfa cold stress.

Figures

Full text

1. Introduction

Hard winters, freezing damage, and recurrent spring cold are challenges frequently encountered by alfalfa (Medicago sativa), and they often occur in high latitudes in the north, constraining forage crop distribution and production [1,2]. Low temperatures cause significant physiological and cellular damage to alfalfa, such as mitochondrial dysfunction, membrane and photosystem injury, and reduced growth. These kinds of damage usually inhibit alfalfa root and shoot activity or growth, carbon metabolism, and photosynthesis [3,4,5,6]. Adaptive mechanisms in alfalfa evolution can reduce their impact.

Since plants are sessile organisms unable to escape from survival crises, they adopt multiple mechanisms to avoid stresses [7]. Plants adapt and survive by readjusting their biochemical activity through gene expression compositions in the face of temperature changes. Genes produce coding RNAs (mRNAs), which function as protein synthesis templates to synthesize or inhibit a substance, influencing how plants sense, respond to, and survive various stresses, from metabolism to precise regulation, and they produce a complex regulatory network [8,9,10,11,12]. In alfalfa’s cold acclimation, to improve its cold tolerance, sucrose is first decomposed then transformed into starch when the photosignal changes, and sucrose is used to synthesis raffinose at the same time; after the temperature cools down, starch eventually hydrolyzes into glucose, and glucose is used to produce sucrose. Therefore, genes such as BAM, SS, SUS, SPS, and RFS were found to be dynamically expressed for stable carbon balance [13]. Research indicates that MsMYB12 in alfalfa can attach directly to the promoter of MsFLS13, which accumulates flavonols and inhibits ROS, enhancing antioxidant capacity under cold stress [14]. In addition, numerous studies have found that ICE-CBF-COR plays an important role in cold stress signaling [15,16]. CBF activates downstream genes COLD REGULATED (COR) and adjusts the physiological and biochemical properties inside plant cells to promote cold tolerance [13,17,18].

Apart from coding RNAs, recent studies have demonstrated that non-coding RNAs (ncRNAs) can affect gene expression, and the competitive endogenous RNA (ceRNA) of these ncRNAs was used to explain the mechanisms underlying the regulation of gene expression [19], influencing plant growth, development, and survival under adverse conditions. Non-coding RNAs (ncRNA), including long non-coding RNAs (lncRNAs), circular RNA (circRNA), and small RNA (sRNA), such as microRNAs (miRNAs), are important members of the RNA family [20]. Recently, evidence suggests that ncRNAs originating from intergenic regions, antisense strands of protein-coding genes, and pseudogenes are involved in regulating transcription, post-transcriptional processes, and mRNA translation [21,22,23,24].

To date, there are few studies on stress-responsive non-coding RNA in alfalfa under cold stress. In this investigation, we discovered mRNAs, miRNAs, lncRNAs, and circRNAs that respond to cold stress in the “Zhaodong” alfalfa variety and built a cold-responsive ceRNA regulatory network. The research provides a theoretical underpinning for future studies on how ncRNAs function in alfalfa under cold stress.

2. Materials and Methods

2.1. Plant Materials and Treatments

The cultivar alfalfa “Zhaodong”, which has excellent cold resistance, was used for whole-transcriptome (WT) analysis and degradome sequencing, with the aim to clarify the role of coding and non-coding RNA in alfalfa cold stress. Zhaodong alfalfa was a local cultivar bred by the Heilongjiang Institute of Animal Husbandry, with excellent cold resistance and 2 fall dormancy level [25,26]. A total of 300 alfalfa accessions (including Zhaodong cultivar) were used for genome-wide association study (GWAS), with the aim to excavate the alfalfa cold resistance gene, then were combined with the results of transcriptome differential gene expression, and candidate regulatory genes were screened.

Seeds of Zhaodong cultivar were first placed on filter papers in sterile Petri dishes at 4 °C in dark conditions, and after germination, they were planted into flowerpots filled with the ratio of 3 peat soil: 1 vermiculite: 1 washed sand, and grown at 25 °C conditions with 16 h in the light and 8 h in the dark cycle. Seeding samples were collected after 3 weeks. In this study, LD refers to the treated group, and YM is the control group, and 10 plant seedlings each from LD and YM were selected randomly as 1 biological replicate, with 3 biological replicates of each group. Seedlings of the LD groups were treated under a 16 light/8 dark photoperiod at −4 °C for 24 h, and seedlings of the YM groups were treated under a 16 light/8 dark photoperiod at 25 °C for 24 h; then, leaf samples were collected and quickly frozen in liquid nitrogen and stored at −80 °C for whole-transcriptome analysis and degradome sequencing.

A total of 300 alfalfa accessions were used for GWAS, the accession information and phenotypic data were showed in Table S1. These accession seeds were collected and provided by the China Academy of Agricultural Sciences Institute of Grassland Research (CAASIGR, Hohhot). Each accession was first bred at the greenhouse; after screening the healthy plants with the same development period, they were transplanted to the field at the 3-leaf stage. Each accession planted 15 plants with 1 m × 1 m distance between the plants, with conventional field management methods. Leaf samples of each accession were collected randomly from 1 plant during the initial flowering period, quickly frozen in liquid nitrogen, and stored at −80 °C. The cold resistance of the accessions was evaluated by the index of regreening time; the earlier the regreening time, the worse the cold resistance of the plant. These phenotypic data were collected an average of continuously 3 years and used for GWAS.

2.2. RNA Extraction and RNA Sequencing

Following the manufacturer’s instructions, total RNA was extracted from each sample using TRIzol reagent (Invitrogen, CA, USA). The quantity and purity of total RNA were assessed using the Bioanalyzer 2100 and the RNA 6000 Nano LabChip Kit (Agilent, CA, USA, 5067-1511), ensuring an RIN number greater than 7.0. Approximately 5 ug of total RNA was used to remove ribosomal RNA, according to the instructions in the Ribo-Zero Gold rRNA Removal Kit (Illumina, San Diego, CA, USA). We carried out the 2 × 150 bp paired-end sequencing (PE150) on an Illumina Novaseq^TM^ 6000 at LC-Bio Technology Co., Ltd., Hangzhou, China, following a quality and integrity assessment and the manufacturer’s protocol. Cutadapt (https://cutadapt.readthedocs.io/en/stable/, version:cutadapt-1.9, accessed on 13 March 2025) was used to further filter the reads obtained from the sequencing machines. Subsequently, we utilized the HISAT2 package (version 2.2.1) to align all sample reads to the Medicago sativa reference genome (https://daehwankimlab.github.io/hisat2/, accessed on 13 March 2025). All sequencing experiments were carried out with 3 biological replicates.

2.3. DNA Extraction, Sequencing, and Genome-Wide Association Study (GWAS)

The DNA from the leaves of 300 alfalfa accessions was extracted using the CWBIO Plant Genome DNA Extraction Kit (Cowin Biosciences, Beijing, China) and sequenced on the BGI-Shenzhen DNBSEQ platform (BGI, Shenzhen, China). Following quality control, the initial sequencing data were matched to the reference genome of the alfalfa cultivar “XinJiangDaYe” [27]. SAMtools software (v1.2) was then used to obtain filtered and sorted BAM files. High-quality SNPs for the GWAS were acquired using GATK Haplotype Caller [28]. These were then applied on the fastlmm model, which can reveal more genetic loci and true positives. The broom R package (v3.2) was utilized to evaluate the significance of differences between genes. The outcomes were displayed as Manhattan and QQ plots generated with the R package qqman (v3.2) [29]. Candidate genes were discovered within a 20 kb area surrounding each significant SNP, as per the reference genome. The identified genes were then annotated in GO and KEGG databases [30,31].

2.4. Detection and Analysis of mRNAs

Differentially expressed (DE) gene or RNA analysis was performed by DESeq^2^ software (v3.2.5) and edgeR (v3.2.2) between the treatment (LD) and control (YM) groups. The fragments per million exons per thousand bases (FRKM) method was used to determine transcript expression levels. The genes or mRNAs with the parameter of false discovery rate (FDR) < 0.05 and absolute fold change ≥ 2 were considered differentially expressed genes or mRNAs [32,33,34]. All DEGs were associated with GO terms in the Gene Ontology database and the KEGG pathway database for annotation and functional enrichment analysis [30,31].

2.5. Identification and Investigation of miRNAs

The TruSeq Small RNA Sample Prep kits (Illumina, San Diego, CA, USA) were used to sequence the sRNA from LD and YM samples, sequencing on the Illumina Hiseq 2000/2500 platform with a single-ended 50 bp (SE50) read length. ACGT101 miR (v4.2, LC Sciences, Houston, TX, USA) was used to process the raw data of miRNA. To acquire clean data, the 3′ adapters and junk reads were first discarded. Subsequently, reads with 18–25 nt base length were aligned to pre-miRNA and the reference genome in miRBase (v22.1) through a BLAST search to detect known miRNAs. The identification of DEmiRNAs was based on normalized deep-sequencing counts and DEGseq, requiring absolute log^2^FC values more than 1 and p-values below 0.05. The GSTAr v1.0 method was used to determine miRNA binding sites and the targeted gene prediction according to the most abundant miRNAs [35,36]. As described in Section 2.4, we performed the functional annotation for the most frequent miRNA targets.

2.6. Detection and Assessment of lncRNAs

The detection and assessment of lncRNAs in LD and YM samples were conducted using the StringTie software (http://ccb.jhu.edu/software/stringtie/, v2.1.6, accessed on 19 March 2025), which was employed to assemble reads for transcript reconstruction. Using gffcompare software (http://ccb.jhu.edu/software/stringtie/gffcompare.shtml, v gffcompare-0.9.8, accessed on 4 April 2025), all reconstructed transcripts were aligned with the reference genome and classified into twelve categories to identify new transcripts. The class codes ‘i, j, o, u, x’ were used to define novel transcripts, which were then employed to predict lncRNA [37]. Initially, transcripts that overlapped with known mRNAs and were shorter than 200 base pairs were removed. Then, CPC v 0.9 (CPC score below 0.5) and CNCI v 2.0 (CNCI score under 0) were used to identify new transcripts with coding potential. The lncRNAs and DElncRNAs expression levels were identified were determined in the same way as DEmRNAs. Python (v2.7) was used to predict the cis-target genes of lncRNAs from coding genes situated within 100,000 bp both upstream and downstream. Next, we conducted GO and KEGG enrichment analysis of lncRNAs target genes, using methods comparable to that outlined in Section 2.4.

2.7. Detection and Examination of circRNAs

In this investigation, CIRCExplorer2 (v 2.2.6) and CIRI (v2.0.2) were used for circRNA identification of LD and YM samples, based on the structural features and splicing characteristics of circRNAs. The unique circRNAs were produced from unmapped reads by tophat-fusion and CIRCExplorer2 or CIRI. The levels of circRNA expression were measured using the split reads per billion mappings (SRPBM) technique, and the DEcircRNA analysis was performed by edgeR software with the parameter of p-value < 0.05 and fold change > 2 [38,39,40]. According to Section 2.4, GO and KEGG enrichment analyses were then conducted on the DEcircRNA host coding genes.

2.8. Identification and Analysis of lncRNA Co-Expression Analysis and ceRNA Network Construction Degradome Sequencing Verification

The prediction of plant ceRNA analysis based on LD and YM samples sequencing is divided into two parts, miRNA-mRNA and miRNA-lncRNA/circRNA. Prediction of circRNAs targets miRNA using TargetScan (5.0) and miranda (3.3a) with TargetScan_score ≥ 50 and miranda_Energy < −10, as circRNAs competitively bind to miRNA and indirectly affect mRNA expression [41,42]. A network visualization was constructed by Cytoscape (v3.7.1) [43]. The constructed library was sequenced using Illumina Hiseq2000/2500. Degradome sequencing was performed by CleaveLand 3 pipeline, using Oligomap short reading frame calibrator finding mRNAs that correspond to the degradation sequence [44,45].

2.9. Statistical Analyses and Figures Creation

Statistical analyses and data were plotted using Origin 2021. Figures and statistical representations were carried out using Adobe Illustrator 2019.

3. Results

3.1. DEmRNAs’ Response to Cold Conditions

A total of six cDNA libraries, including those from “Zhaodong” alfalfa leaves of YM samples and LD samples, generating more than 77,629,034 raw reads through RNA sequencing and at least 75,936,528 clean reads, were preserved for further examination (Table S2). The reads were aligned to the alfalfa genome, achieving an average mapping ratio of over 80.15%. Q20 values were 99.98%, and Q30 values exceeded 98.47%. The GC content was more than 40%, indicating high-quality reads that could be used for differential gene expression analysis. More than 109,110 genes and 132,503 transcripts were successfully annotated in the six databases. Furthermore, we conducted the correlation analysis between samples of these libraries based on expression levels. The biological repeatability between samples is good according to the combination of the correlation heat map and the PCA map (Figure S1).

In total, 12,693 mRNAs were identified as expressed genes according to FPKM (Table S3). The volcano plot revealed 46,936 differential expression mRNAs (DEmRNAs) in the LD group (Figure 1A), while a total of 8963 were up-regulated and 3730 were down-regulated in response to cold stress (Table S4). GO annotation analysis identified the cellular component linked to alfalfa’s response to cold stress, including chloroplast (GO:0009507, 1736 DEmRNAs), cytosol (GO:0005829, 1245 DEmRNAs), mitochondrion (GO:0005739, 1162 DEmRNAs), chloroplast stroma (GO:0009570, 488 DEmRNAs), and cell wall (GO:0005618, 437 DEmRNAs) (Figure 1B, Table S5). The chloroplast could synthesize carbohydrates, and carbon catabolism and generation of energy took place in the mitochondria. Regarding biological processes, DEmRNAs were predominantly associated with the oxidation–reduction process (GO:0055114, 476 DEmRNAs), response to cold (GO:0009409, 273 DEmRNAs), response to salt stress (GO:0009651, 267 DEmRNAs), response to cadmium ions (GO:0046686, 245 DEmRNAs), and translation (GO:0006412, 218 DEmRNAs), but they were most significantly enriched in the starch catabolic process (GO:0005983, 40 DEmRNAs, p-value = 9.88 × 10^−24^), which included 36 up-regulated and 4 down-regulated DEmRNAs. Regarding molecular function, DEmRNAs were associated primarily with structural constituents of ribosomes (GO:0003735, 290 DEmRNAs), copper ion binding (GO:0005507, 128 DEmRNAs), peroxidase activity (GO:0004601, 110 DEmRNAs), heme binding (GO:0020037, 109 DEmRNAs), and unfolded protein binding (GO:0051082, 83 DEmRNAs).

Significant enrichment of DEmRNAs in pathways was indicated by the KEGG enrichment analysis, including ribosome (ko03010, 326 DEmRNAs, p-value = 0.00 × 10^0^); amino sugar and nucleotide sugar metabolism (ko00520, 184 DEmRNAs, p-value = 5.81 × 10^−10^); ribosome biogenesis in eukaryotes (ko03008, 143 DEmRNAs, p-value = 1.38 × 10^−9^); circadian rhythm–plant (ko00270, 149 DEmRNAs, p-value = 1.45 × 10^−7^); alanine, aspartate, and glutamate metabolism (ko00250, 62 DEmRNAs, p-value = 2.34 × 10^−6^); and starch and sucrose metabolism (ko00500, 241 DEmRNAs, p-value = 7.07 × 10^−6^) (Figure 1C, Table S6). Among these KEGG terms, ribosome and starch and sucrose metabolism pathway enriched the most number of DEmRNAs. Taken together, these results show that in the reprogramming of critical molecular pathways, DEmRNAs involved in cold resistance are focused mostly on systematic adjustment of translational capacity and precise redistribution of carbon metabolism.

Through enriched GO and enriched KEGG functions, we found that alfalfa responses to cold stress can be summarized to sustain cellular energy supply, and carbon reallocation (starch breakdown) was central to the acute energy strategy. A total of 194 DEmRNAs associated with the starch and sucrose metabolism pathway were up-regulated and 47 DEmRNAs were down-regulated (Table S7). In this term, two AMY3 (MS.gene63187 and MS.gene062524), two EG6 (MS.gene069582 and MS.gene55793), and one PGSIP8 (MS.gene71946) were the top five significantly enriched DEgenes, whereas MS.gene003205 (EG25), MS.gene56743 (GPI), MS.gene34751 (AMY3), MS.gene09696 (BAM3), and MS.gene015359 (BAM3) were the top five highly expressed up-regulated DEgenes. In addition, there is a multi-level collaborative relationship between pathways of amino sugar and nucleotide sugar metabolism and starch and sucrose metabolism; they shared the same intermediate materials or exchanged materials (such as UDP-glucose), and both of them work in plant response against abiotic stress. As for the pathway of amino sugar and nucleotide sugar metabolism, with 128 up-regulated DEmRNAs and 32 down-regulated DEmRNAs (Table S7), 3 UAM1 (MS.gene53818, MS.gene40842, and MS.gene44706) and 2 BXL (MS.gene86920 and MS.gene35614) were the top 5 significantly enriched DEgenes, whereas MS.gene56743 (G6PI1), MS.gene015702 (PGM), MS.gene051339 (BXLAF2), MS.gene62173 (GMP), and MS.gene65895 (BXL7) were the top 5 highly expressed up-regulated DEgenes. These functional DEgenes were the potential candidate genes under cold stress.

Additionally, we analyzed several important transcription factor (TF) families (Figure 1D; Table S8), including the NAC, MYB, WRKY, bHLH, AP2/ERF, MADS-box, zinc finger, GRAS, ARF, and TCP. The zinc finger, MYB, bHLH, and WRKY families comprised the largest groups, with 127, 59, 39, and 21 members, respectively, under cold stress conditions. The zinc finger family member contains at least 20 different zinc finger domains. Among DEmRNAs of these genes, 78 zinc finger family members, 27 MYB family members, 21 bHLH family members, and 17 WRKY family members were up-regulated. A total of ten DEmRNAs of zinc finger were annotated to the indole alkaloid biosynthesis pathway (K00901), and nine were annotated to circadian rhythm–plant (K04712). Four DEmRNAs of MYB were annotated to plant hormone signal transduction (K04075); ten DEmRNAs of bHLH were annotated to the MAPK signaling pathway–plant (K04016) and plant hormone signal transduction (K04075); and nine DEmRNAs of WRKYs were annotated to MAPK signaling pathway–plant (K04016) and plant–pathogen interaction (K04626), which were the top gene annotation numbers of the KEGG pathway (Table S8). It is indicated that the zinc finger family possibly participates in alfalfa cold resistance through plant biological clock regulation as well as secondary metabolite production for oxidation-resistant and cell stability maintenance. The MYB family and bHLH family possibly participate in plant hormone signal transduction to regulate expression of cold-resistant genes. In addition, the bHLH family and WRKY family possibly participate in MAPK signaling to receive the upstream cold signal, then phosphorylate the critical cold-resistant transcription factors, and finally activate the cold tolerance gene expression. The WRKY family might also participate in the protection of bacterial invasion caused by freezing injury.

3.2. GWAS Identifies Significant SNPs Associated with Cold Stress Resistance

To obtain the candidate genes related to cold stress resistance, 300 alfalfa accessions were used for genome-wide association studies (GWASs). Germplasm information is showed in Table S1. To obtain more reliable SNPs, we used GATK (Genome Analysis Toolkit) to identify potential SNP loci, and MLM methods were used in this study. After screening with the threshold p-value < 2.91 × 10^−7^ to identify significant association signals, four significant SNP loci were detected for cold stress resistance (Figure 2). A total of 118 putative genes near the 4 associated loci were identified, and they are presented with gene annotations including possible functions in Table S9. Integrating GWAS and transcriptome analysis to identify the overlapped DEGs, the only one candidate gene with differential expression was MS.gene53818 (Figure 2; Table S9). On the basis of the above results, MS.gene53818 (MsUAM1) was found to be one of significantly enriched DEgenes (Table S4).

3.3. Response of DEmiRNAs to Cold Stress

MicroRNAs are crucial components in regulatory networks, influencing both transcriptional and post-transcriptional processes. The small RNA-seq library obtained more than 10,233,316 clean reads (Table S2). A total of 1965 miRNAs were expressed in alfalfa (Table S10); these could be classified into five groups, including 481 known miRNAs (gp1), 507 conservative miRNAs (gp2a, gp2b), and 926 potential novel miRNAs (gp4). The identified miRNAs typically ranged from 18 to 25 nt in length, and 46.62% had a length of 24 nt and 30.70% a length of 21 nt (Figure 3A). The known miRNAs were categorized into 86 families, among which 58 families contain at least 3 members and 14 families contain only 1 member. The top 10 families with the highest number of members were MIR156 (43 miRNAs), MIR159 (35 miRNAs), MIR166 (30 miRNAs), MIR2655 (29 miRNAs), MIR167_1 (24 miRNAs), MIR2592 (24 miRNAs), MIR169_2 (20 miRNAs), MIR171_1 (16 miRNAs), and MIR482 (13 miRNAs) (Figure 3B, Table S11).

A total of 149 and 91 mRNAs were specific to the LD and YM groups, respectively (Figure 3C). A total of 223 DEmiRNAs were further identified (p-value <0.01). Compared with those in the control group YM, 134 were up-regulated and 89 down-regulated in the LD group. Of these, 39 DEmiRNAs were highly expressed. Under cold stress, the miRNAs PC-3p-177_6262, gma-miR4416c-5p_R+1_1ss2TC (MIR4416), gma-miR319f, PC-5p-982 1440, and mtr-miR396b-3p (MIR396) showed the highest levels of expression, whereas PC-5p-76_16535, gma-MlR5770b-p5_1ss1TA, stu-MIR6026-p3_2ss7TA19CA (MIR6026), and mtr-MIR2671h-p3 (MIR2655) presented the lowest expression levels under cold stress (Figure 3D, Table S12).

A total of 2524 putative mRNAs targeted by DE-miRNAs were predicted (Table S13). Almost all miRNAs, except for 11 miRNAs, targeted more than 1 coding gene, and 12 miRNAs targeted more than 50 coding genes. To further investigate the functions of the cold-responsive DEmiRNAs, GO and KEGG functional enrichment analyses were conducted on the DEmiRNAs target genes. The results revealed that these target genes were enriched in 696 GO terms (Table S14). Of the 137 enriched GO terms of the cellular component category, the top three significantly enriched terms were nucleus (GO:0005634), SCF ubiquitin ligase complex (GO:0019005), and catalytic step 1 spliceosome (GO:0071012) (Figure 4A, Table S14). A total of 514 GO categories were enriched in the biological process category (Figure 4B, Table S14); terms of DNA integration (GO:0015074), anther development (GO:0048653), and leaf development (GO:0048366) were the top three significantly enriched terms. Of the 318 enriched GO terms in the molecular function category (Figure 4C, Table S14), the terms of RNA-directed DNA polymerase activity (GO:0003964), DNA-directed DNA polymerase activity (GO:0003887), and DNA-binding transcription factor activity (GO:0003700) were significantly enriched. KEGG enrichment analysis revealed that the target DEmRNAs of the DEmiRNAs were mostly annotated in plant hormone signal transduction (62 target DEmRNAs) (Figure 4D, Table S15), protein processing in endoplasmic reticulum (52 target DEmRNAs), and starch and sucrose metabolism (37 target DEmRNAs), but indole alkaloid biosynthesis was the most significantly enriched pathway with highest p-value (0.000002), followed by plant hormone signal transduction (p-value = 0.00004), and starch and sucrose metabolism ranked 11th (Figure 4E, Table S15).

3.4. Response of DElncRNAs to Cold Stress

lncRNAs usually act as cis-regulators of nearby genes and regulate their transcription and expression. To further understand the regulatory roles of alfalfa lncRNAs in response to cold stress, we analyzed the identified lncRNAs, DElncRNAs, and their target genes. A total of 10,751 lncRNAs were identified from alfalfa, including 1852 significant expressed lncRNAs with 819 up-regulated and 1033 down-regulated (Figure 5A). Of these 1523 DElncRNAs, across all samples, intergenic lncRNAs (u-type) were the most common, with 141 intronic lncRNAs (i-type), 65 antisense lncRNAs (x-type), 64 sense lncRNAs (o-type), and 59 bidirectional lncRNAs (j-type) (Table S16), which suggest that lncRNAs are mainly derived from the non-coding regions. We further compared the structural features of lncRNAs and found that the majority of lncRNAs were shorter than 900 nt, and most of these lncRNAs contained one or two exons (1726) (Figure S2). The volcano plot highlights that the expression levels of the lncRNAs MSTRG.72.1 (u-type), MSTRG.35733.1 (u-type), MSTRG.79.1 (u-type) and MSTRG.34811.1 (i-type) were significantly up-regulated, whereas the expression levels of the lncRNAs MSTRG.66445.1 (u-type), MSTRG.44038.1 (u-type), MSTRG.54658.1 (u-type), MSTRG.458.1 (u-type), MSTRG.7081.1 (u-type), and MSTRG.66848.1 (u-type) were significantly down-regulated under cold stress (Figure 5A, Table S16).

The 542 DElncRNAs were predicted to have potential cis-regulatory effects on 695 DEmRNAs. A total of 405 DEmRNAs and DElncRNAs presented one-to-many matches, whereas 379 DEmRNAs and DElncRNAs presented one-to-one matches (Table S17). The findings revealed that 4350 target genes were enriched in 2839 GO terms. Central vacuole (GO:0042807) and plant-type vacuole membrane (GO:0009705) were the most enriched cellular component terms, while response to auxin (GO:0009733) and ion transmembrane transport (GO:0034220) were the most enriched biological process terms; water channel activity (GO:0015250) and glycerol channel activity (GO:0015254) were the most enriched molecular function terms (Figure 5B, Table S18). KEGG enrichment results demonstrated that 1752 target DEmRNAs of DElncRNAs were enriched in 134 metabolic pathways; among them, ascorbate and aldarate metabolism (ko00053, p-value = 0.0009, 39 DEmRNAs), plant hormone signal transduction (ko04075, p-value = 0.0019, 157 DEmRNAs), flavone and flavonol biosynthesis (ko00944, p-value = 0.0028, 22 DEmRNAs), and glycosylphosphatidylinositol (GPI)-anchor biosynthesis (ko00563, p-value = 0.0036, 16 DEmRNAs) were the most enriched pathways (Figure 5C, Table S19). These data suggest that under cold stress, DElncRNAs and DEmRNAs may play a key role in the regulation of plant growth and development during environmental adaptation.

3.5. Response of DEcircRNAs to Cold Stress

CircRNAs usually act in cis/trans patterns on target genes and frequently compete with mRNAs for miRNA binding, thereby modulating mRNA degradation and translation. Employing two methods for subsequent analyses (Table S20), we identified—along with the mRNAs, miRNAs, and lncRNAs—258 circRNAs, including 129 exonic circRNAs (circRNA type), 54 intronic circRNAs (ciRNA type), and 75 intergenic circRNAs (intergenic type) (Figure 6A), which suggests that circRNAs are mainly derived from the coding regions of genes. Most genes have only one circRNA, and only one gene (MS.gene057745) harbors more than one circRNA (circRNA153, circRNA154) (Figure 6B, Table S20). Almost all the circRNAs are distributed across all chromosomes; chromosome 3.2 contains the most, while chromosome 2.2 contains only one circRNA (Figure 6C, Table S20). The vast majority of circRNAs contain one to two exons, and five circRNAs have over ten exons (Figure 6D). Among all circRNAs, thirteen DEcircRNAs were identified in alfalfa’s response to cold stress (Figure 6E, Table S21), with six up-regulated and seven down-regulated. The heat map of the DEcircRNAs was generated to display the DEcircRNA expression profiles in the individually treated samples (Figure S3). A total of 11 DEcircRNAs were predicted to exert effects on 74 DEmiRNAs (Table S22). The results show that circRNA27, circRNA33, circRNA83, circRNA38, circRNA75, circRNA82, and circRNA106 exhibited over three connections, which is more than other circRNAs did, pointing to their potential importance as ceRNAs within the network.

3.6. CeRNA Network Regulation in Alfalfa Reaction to Cold Stress

To clarify the worldwide regulatory network of ncRNAs and mRNAs associated with cold stress in “Zhaodong” alfalfa, a ceRNA network of key DEmRNAs involved in alfalfa’s resistance to cold stress was established. The ceRNA network comprised 98 DEmRNAs, 274 DElncRNA, 33 DEmiRNA, and 8 DEcircRNA (Tables S23 and S24). These DEgenes were significantly enriched in starch and sucrose metabolism, amino sugar and nucleotide sugar, and plant hormone signal transduction KEGG pathways. Within the ceRNA network, 23 DEmRNAs, 11 DEmiRNAs, 28 DElncRNAs, and 8 DEcircRNAs were validated by degradome sequencing (Figure 7).

The ceRNA network comprised several proteins and transcription factors that play an important role in alfalfa resistance to cold stress, including the calmodulin (CaM), rubisco accumulation factor 1.1 (RAF1.1), berberine bridge enzyme-like 13 (BBE-like 13), ARM repeat superfamily protein (ARM), CROWDED NUCLEI 1 isoform X1 (CRWN1), branched-chain amino acid aminotransferase (BCAT), glucan endo-1,3-beta-glucosidase (BGLU), fumarate hydratase (FH), cysteine proteinase (COT44), and transcription factor WRKY33 (Tables S23–S25). The findings revealed that three CaMs (MS.gene016892, MS.gene47025, MS.gene37614) and one RAF1.1 (MS.gene74997) were regulated by stu-MIR6026-p3_2ss7TA19CA, MSTRG.51140.1 and MSTRG.28735.1. BBE-like 13 (MS.gene35113) and ARM (MS.gene065493) were regulated by mtr-MIR2592ab-p5. CRWN1 (MS.gene43202) and WRKY33 (MS.gene049457, MS.gene41592) were regulated by PC-3p-79740_34. BCAT (MS.gene77492) was regulated by vvi-miR845a_L-1R-2_1ss15AC. FHs (MS.gene20767, MS.gene07245) were regulated by PC-3p-58625_51. COT44 was regulated by mtr-miR396a-5p. The results indicate that the ceRNA network may enhance the cold stress resistance of alfalfa from multiple dimensions.

4. Discussion

The development and growth of alfalfa are significantly influenced by cold stress. Nowadays, many studies have delved into the roles and regulatory processes of genes involved in cold stress. Even though numerous protein-coding genes associated with cold stress have been functionally characterized, the regulatory roles of non-coding RNAs in alfalfa facing cold stress are still not well understood. In this investigation, we discovered cold-responsive DEmRNAs, DEmiRNAs, DElncRNAs, and DEcircRNAs in alfalfa leaves under cold stress and constructed a ceRNA network responsive to cold, focusing on potential RNA interactions.

4.1. Carbon Balance Reconstruction Is the Key Strategy for Alfalfa in Cold Stress and UAM1 Was Identified as the Critical Candidate Gene

Many alfalfa cultivars are highly vulnerable to abiotic stressors in the early development or regeneration period. The abrupt temperature change is a key factor limiting growth or plant regeneration in alfalfa seedings. This study comprehensively demonstrated the complex transcriptomic reprogramming of alfalfa under cold stress through high-quality transcriptome sequencing. Our study not only identified many differentially expressed genes (DEGs) but also converted the gene list into biologically meaningful pathways and functional modules by in-depth GO and KEGG enrichment analyses. This helps clearly delineate the core strategies of alfalfa in responding to cold stress: carbon metabolism redistribution and energy supply adjustment.

The starch and sucrose metabolism pathway is the center of plant carbon metabolism in alfalfa response to cold stress [13,46]. Starch and sucrose metabolism was also found to be one of the most significantly enriched terms (Figure 1C, Table S6), with the highest number of DEmRNAs, and it is mainly regulated by up-regulated DEGs, which indicates that the conversion between starch and sucrose (starch hydrolyzed to maltose then into glucose, followed glucose was used to synthesis sucrose) is an early-stage core strategy for alfalfa to cope with cold stress. We noticed that several significant DEGs were involved in starch hydrolysis, decomposition, and reconstruction of the cell wall. The rapid degradation of starch into monosaccharides provides important carbon skeletons, signaling substances, osmoprotectants for other metabolic pathways, and substrates for plant cell wall remodeling [47,48,49,50,51].

Moreover, by integration of GWAS and mRNA analysis, we identified an important candidate gene associated with alfalfa’s cold tolerance, UDP-arabinopyranose mutase 1 (MsUAM1) (Figure 2, Table S9). MsUAM1 was up-regulated in “Zhaodong” alfalfa cold resistance (Tables S7 and S9). Arabidopsis mur5 is an allele of RGP2, which functions on interconverts UDP-Arap and UDP-Araf, and either the missense or null mutant of mutant mur5 exhibited UAM activity deficiency and the low cell wall UDP-Ara phenotype [52]. UDP-arabinopyranose mutase is essential for plant cell wall biosynthesis and development. UAM catalyzes the reversible interconversion of UDP-Arap (pyranose form) and UDP-Araf (furanose form), with the furanose form being the biologically active donor for arabinose incorporation into cell wall polymers, which connects the metabolism of primary carbohydrates (starch and sucrose) to the biosynthesis of cell wall polysaccharides [53,54,55,56]. Starch and sucrose metabolism generate UDP-glucose (UDP-Glc), a central nucleotide sugar. UDP-Glc is a precursor for several other UDP-sugars and acts through a series of interconversion reactions [51,57]. UAMs are required for the biosynthesis of arabinan side chains in pectins and arabinoxyloglucans, which are crucial for cell wall integrity, leaf expansion, and stress responses [53,54,55,56]. Cell wall arabinose content and structure, regulated by UAM, may influence plant tolerance to cold. Proper arabinose incorporation into cell wall polymers can affect wall flexibility and porosity, which are important for withstanding the physical stresses of freezing and thawing [53,55], although direct evidence for specific action of UAMs in cold stress is still limited.

4.2. Cold-Responsive ncRNAs in Alfalfa

Whole-transcriptome analysis enables comprehensive profiling of gene expression in crops under stress, revealing the complex regulatory networks and key genes involved in stress adaptation and resilience. Based on our sequencing results, 134 DEmiRNAs were up-regulated and 89 DEmiRNAs were down-regulated by cold stress in alfalfa (Figure 3C). The vast majority of miRNAs targeted more than one coding gene (Table S13), and this result reflects the complexity of their regulatory network. Previous research supports the finding that most plant miRNAs target multiple coding genes, enabling broad and coordinated regulation of plant gene expression. In Triticum aestivum, 524 targets were identified for 124 miRNA families, averaging over 4 targets per miRNA, with some miRNAs targeting even more genes. A high degree of target multiplicity was observed in Panicum virgatum, where 121 miRNAs were predicted to target 839 protein-coding genes, similarly to Pseudostellaria heterophylla [58,59,60]. This multi-targeting enables miRNAs to coordinate the regulation of gene families involved in key processes like development, stress response, and metabolism [61,62]. Indeed, our results show that targets of DEmiRNAs were enriched in various KEGG pathways such as plant hormone signal transduction, protein processing in endoplasmic reticulum, and starch and sucrose metabolism (Figure 4E, Table S15).

We identified 1852 DElncRNAs and predicted their target genes, finding that these DElncRNAs may play a crucial role in the adaptation of alfalfa to cold stress by regulating key genes related to cell homeostasis, hormone signaling, antioxidant metabolism, and cell wall modification, as indicated by GO terms such as central vacuole (GO:0042807), plant-type vacuole membrane (GO:0009705), water channel activity (GO:0015250), and response to auxin (GO:0009733) (Figure 5B, Table S18). In plant cells, vacuoles are the predominant organelles, frequently comprising up to 90% of the cell’s volume, and their primary roles involve the following: turgor pressure maintenance for supporting plant structure and development; storage of metabolites for supporting metabolism; detoxification and defense against pathogens; and cellular homeostasis and stress response to compartmentalize harmful substances and degrade macromolecules [63,64,65,66]. Low temperatures could result in cellular dehydration as well as larger ice crystals. Aquaporins play a crucial role in plant cold stress tolerance by maintaining water balance, enhancing antioxidant defenses, and regulating stress-responsive genes. The up-regulation of aquaporin genes increases hydraulic conductivity and water transport capacity, which is essential for low-temperature environment adaptation [67,68,69]. Notably, target genes were also enriched in ascorbate and aldarate metabolism (ko00053) which is the most enriched pathway (Figure 5C, Table S19). Enhancing ascorbic acid biosynthesis with oxidative damage reduction improved alfalfa’s tolerance to not only cold but also saline-alkali stress [70]. GO and KEGG enrichment of the target genes of DElncRNAs slightly differed from DEmiRNAs’ target genes, with more enrichment in those specific metabolic enzymes or channel proteins with functional execution. It seems that lncRNAs could function not only in transcription and post-transcriptional processes but also in direct gene regulators [23,71,72].

The ceRNA hypothesis proposes that mRNAs, lncRNAs, and circRNAs can regulate each other by competing for shared miRNAs and that the competition forms complex regulatory networks that add a new dimension to post-transcriptional gene regulation. A total of 413 nodes (98 DEmRNAs, 274 DElncRNAs, 33 DEmiRNAs, 8 DEcircRNAs) were identified, but only 70 of them (23 DEmRNAs, 11 DEmiRNAs, 28 DElncRNAs, 8 DEcircRNAs) were validated by degradation sequencing (Tables S23 and S24, Figure 7). The genes in the network were significantly enriched in plant hormone signal transduction, plant–pathogen interaction, starch and sucrose metabolism, and ascorbate and aldarate metabolism; this corroborates the preceding analysis of DEmRNA, DEmiRNA, and DElncRNA, which suggested that these are the core strategies of alfalfa’s response and adaptation to cold stress, finely controlled by complex ceRNA regulatory networks. In plants, miRNAs can inhibit translation initiation by binding to complementary sites on target mRNAs, but binding can sterically block ribosome recruitment or movement, repressing protein synthesis without necessarily degrading the mRNA [73,74,75]. In this ceRNA, we predicted most ncRNAs may act as miRNA sponges, but it is worth noting that PC-3p-79740_34 regulates the expression of two WRKY33 (MS.gene049457, MS.gene41592) genes, and two lncRNAs (MSTRG.48783.1, MSTRG.14777.1) binding to PC-3p-79740_34 might act as miRNA sponges and reduce the inhibition of WRKY33. WRKY33 has been reported to be associated with cold in Solanum lycopersicum, Vitis amurensis, and Paeonia lactiflora [76,77,78]. Overall, these results indicate that ceRNA networks are involved in regulating gene expression to adjust carbon balance, maintain cell homeostasis, and regulate metabolism. Thus, the findings enrich our understanding of how ceRNAs might regulate alfalfa cold tolerance and provide genetic references for further analysis of cold stress response mechanisms.

5. Conclusions

We applied WTRNA-seq to explore the response to cold stress in alfalfa leaves, identifying 46,936 DEmRNAs, 223 DEmiRNAs, 1852 DElncRNAs, and 13 DEcircRNA. Our findings indicate that specific lncRNAs and circRNAs serve as ceRNAs in alfalfa’s response to cold conditions. A ceRNA network comprising 28 DElncRNAs, 8 DEcircRNAs, 11 DEmiRNAs, and 23 DEmRNAS triplets was constructed. The involvement of competitive lncRNA/circRNA/miRNA-mRNA regulatory networks is seen in amino sugar and nucleotide sugar metabolism, starch and sucrose metabolism, and plant hormone signal transduction. Regulating carbon balance is a crucial method for alfalfa to enhance its cold resistance. WTRNA-seq and GWAS identified several candidate genes, including UAM1 and WRKY33, that play a role in the cold stress response. In summary, these findings establish a basis for examining ncRNA reactions to cold stress in alfalfa breeding.