Widely Targeted Metabolomics and Transcriptomics Analysis of the Response and Adaptation Mechanisms of Trifolium ambiguum to Low-Temperature Stress

Abstract

Caucasian clover (Trifolium ambiguum M.Bieb.) is a perennial legume known for its exceptional cold tolerance, commonly used in agriculture and ecosystems in cold climates. Given the impact of climate change, enhancing the cold adaptation of Caucasian clover is crucial for sustaining agricultural productivity. This study employs metabolomics, transcriptomics, and Weighted Gene Co-expression Network Analysis (WGCNA) to investigate the molecular mechanisms of Caucasian clover’s response to low-temperature stress. Metabolomic analysis showed that low-temperature stress triggered the accumulation of fatty acids, amino acids, and antioxidants, which are critical for maintaining membrane stability and antioxidant capacity, thus protecting the plant from oxidative damage. Transcriptomic analysis revealed significant upregulation of genes involved in cold adaptation, particularly those related to antioxidant defense, membrane lipid repair, and signal transduction, including genes in the ABA signaling pathway and antioxidant enzymes, thereby improving cold tolerance. WGCNA identified gene modules closely linked to cold adaptation, especially those involved in antioxidant defense, fatty acid metabolism, signal transduction, and membrane repair. These modules function synergistically, with coordinated gene expression enhancing cold resistance. This study also investigated the isoflavonoid biosynthesis pathway under low-temperature stress, highlighting its role in enhancing antioxidant capacity and cold tolerance. Low-temperature stress induced upregulation of key enzyme genes, such as Isoflavone Synthase (IFS) and Isoflavone-7-O-Glucosyltransferase (IF7GT), promoting antioxidant metabolite accumulation and further enhancing the plant’s cold adaptation. Overall, this study offers novel molecular insights into the cold tolerance mechanisms of Caucasian clover and provides valuable theoretical support for breeding cold-resistant crops in cold climates.

1. Introduction

Caucasian clover (Trifolium ambiguum M.Bieb.), a perennial legume, demonstrates remarkable environmental adaptability, primarily due to its enhanced cold resistance [1,2]. Compared to other legumes such as alfalfa and white clover, Caucasian clover is uniquely characterized by its superior cold tolerance, high nutritional value, and substantial protein content, making it an essential forage resource for cold regions [3]. Its adaptability and resilience under extreme environmental conditions highlight its potential for sustainable agricultural practices and livestock production. Furthermore, the economic benefits of cultivating Caucasian clover are significant, as it contributes to increased forage yield and quality in cold regions, supporting both livestock and agricultural sustainability.

Low-temperature stress is a significant factor limiting plant growth and productivity [4,5]. In low-temperature environments, the stability of plant cell membranes, metabolic activity, and physiological functions are profoundly affected, resulting in oxidative damage, metabolic disturbances, and reduced photosynthetic efficiency [6]. These alterations directly impact plant growth, development, and stress tolerance, particularly under cold conditions in early spring and autumn, posing substantial challenges to agricultural and pastoral productivity [7,8,9]. Therefore, understanding the molecular mechanisms underlying plant responses to low-temperature stress not only illuminates their adaptive strategies but also paves the way for enhancing crop cold tolerance through molecular breeding, thereby creating new opportunities for agricultural production in colder regions.

In recent years, advancements in omics technologies have facilitated the widespread application of transcriptomics and metabolomics in research on plant stress responses [10,11,12]. These approaches have proven invaluable for understanding the molecular basis of crop adaptation to extreme conditions, including cold stress, as demonstrated in recent studies on legumes and model crops such as alfalfa and white clover. Transcriptomic analysis elucidates the dynamic alterations in gene expression and regulatory networks of plants subjected to stress conditions [13], while metabolomics provides deeper insights into how plants respond to environmental changes at the cellular metabolic level [14]. The integrated analysis of these two omics approaches offers a comprehensive view of the multi-layered regulatory mechanisms that plants employ to manage stress. For instance, the combined transcriptomic and metabolomic analysis of Nicotiana tabacum under cold stress revealed the regulatory patterns of lignin and polyphenol biosynthesis pathways, underscoring their significant role in enhancing plant cold tolerance [15]. Additionally, in the study of low-temperature adaptation mechanisms in Euonymus japonicus, alterations in extracellular matrix metabolites under cold conditions were found to enhance the plant’s frost tolerance [16]. Research on different Malus pumila varieties reveals that cold-tolerant ones respond to low temperatures by activating specific metabolic pathways, including those related to carbohydrate and amino acid metabolism, as well as the biosynthesis of phenolic compounds. These pathways may play a pivotal role in enhancing cold resistance [17]. Furthermore, studies suggest that abscisic acid (ABA) signaling and carbohydrate metabolism pathways are essential for cold defense in Argyranthemum frutescens. The plant improves cold tolerance by upregulating relevant genes and accumulating sugars through ABA signaling. Comprehensive transcriptomic and metabolomic analyses of wheat have also revealed the significant influence of low-temperature stress on the abscisic acid (ABA) and jasmonic acid (JA) signaling pathways. Additionally, the proline biosynthesis pathway is critical for low-temperature adaptation, reflecting the regulatory networks involving various cold stress-related genes and metabolites [18].

Although the metabolic and gene expression patterns of various plants under low-temperature stress have been thoroughly investigated, the response mechanisms of Caucasian clover to low-temperature stress remain to be systematically elucidated. The molecular mechanisms underlying its low-temperature tolerance are still in the exploratory phase, with significant gaps in the integration of multi-omics analyses. Therefore, a comprehensive investigation into the low-temperature stress response mechanisms of Caucasian clover will not only provide valuable insights for enhancing its cold tolerance but also offer a reference framework for breeding other cold-tolerant crops.

This study hypothesizes that the response mechanisms of Caucasian clover to low-temperature stress are mediated by specific metabolites and regulatory genes, which play critical roles in enhancing cold tolerance. To test this hypothesis, we employed an integrated approach combining widely targeted metabolomics and transcriptomics, aiming to identify key metabolites and genes involved in the cold stress response. The experimental plan includes analyzing the dynamics of metabolite accumulation and gene expression under low-temperature stress, with particular focus on antioxidant defense, fatty acid metabolism, and membrane repair. By identifying and characterizing these molecular players, this study will provide valuable insights into the genetic and biochemical pathways governing cold tolerance and offer potential candidate genes and molecular markers for breeding cold-resistant varieties. The findings of this study have broad implications, not only advancing the scientific understanding of plant cold tolerance but also contributing to practical strategies for improving agricultural productivity in cold climates.

2. Materials and Methods

2.1. Experimental Materials and Seed Pre-Treatment

This study utilized Caucasian clover seeds as the experimental material, supplied by Inner Mongolia Agricultural University (registration number: N010). To ensure uniformity in germination and the reproducibility of the experiment, the seeds underwent the following pre-treatment steps: (1) Seed cleaning: The Caucasian clover seeds were surface-sterilized with a 0.1% sodium hypochlorite solution for 5 min to eliminate surface impurities and prevent contamination. This was followed by three washes with sterile distilled water to ensure the seeds were free from contamination. (2) Seed washing: The Caucasian clover seeds were rinsed with distilled water to eliminate surface impurities, ensuring that the seeds were free from contamination. (3) Soaking treatment: The cleaned seeds were wrapped in moist filter paper, placed in a Petri dish, and kept moist. They were subjected to a 12 h soaking treatment at 25 °C in the dark to promote uniform germination. (4) Germination culture: The soaked seeds were transferred to a medium containing a 1/2 Hoagland nutrient solution, placed in a Percival incubator, and germinated at 4 °C for seven days under a photoperiod of 16 h of light (light intensity 70%) and 8 h of darkness. (5) Seedling transplanting: After seven days of germination, the seedlings were transplanted into 10 cm × 10 cm culture containers, with the substrate consisting of a 1:1 mixture of vermiculite and Danish Pindstrup soil. The seedlings were cultured for an additional seven days under the same light conditions and environmental controls to facilitate further growth.

2.2. Low-Temperature Stress Treatment

Fourteen-day-old Caucasian clover seedlings were divided into treatment and control groups: (1) Treatment group: Seedlings were subjected to low-temperature treatment at 4 °C and sampled after 2, 6, and 12 h. (2) Control group: Seedlings were maintained at 25 °C without low-temperature treatment. Each group was replicated three times biologically. Immediately following treatment, leaf samples were collected, frozen in liquid nitrogen, and stored at −80 °C for further analysis. Samples were labeled as follows: control group (A), 2 h low-temperature treatment (B), 6 h low-temperature treatment (C), and 12 h low-temperature treatment (D).

The selection of sampling times (2, 6, and 12 h) is based on previous studies of cold stress responses, which used similar time points and considered variations across different stress stages (e.g., early response and adaptation). Preliminary experimental results also confirmed that these time points effectively capture physiological and molecular changes in plants under cold stress, justifying their selection for sampling.

2.3. Widely Targeted Metabolomics Analysis and Data Processing Methods

Widely targeted metabolomics analysis was performed using a UPLC-ESI-MS/MS system (UPLC: Waters Acquity I-Class PLUS; MS: Applied Biosystems QTRAP 6500+), employing a Waters HSS-T3 column (1.8 µm, 2.1 mm × 100 mm) with a gradient program of solvent A (pure water with 0.1% formic acid and 5 mM ammonium acetate) and solvent B (acetonitrile with 0.1% formic acid). The MS parameters included a source temperature of 550 °C, ion spray voltage of 5500 V (positive) and −4500 V (negative), and optimized DP and CE for individual MRM transitions, with nitrogen as the collision gas. After normalizing the original peak areas, principal component analysis (PCA) and Spearman correlation analysis were used to evaluate sample repeatability and quality control. Identified compounds were annotated using KEGG, HMDB, and Lipidmaps databases. Differential metabolites were screened using fold change (FC > 1), T-test p-value (<0.05), and VIP value (>1) from OPLS-DA modeling with 200 permutation tests to validate the model, and KEGG pathway enrichment significance was analyzed using a hypergeometric distribution test. This comprehensive workflow ensures the accurate detection, reliable modeling, and robust identification of differential metabolites.

Metabolomic data were standardized using the internal standardization method, normalizing the total peak area across all samples. This method minimized biases between experiments, ensuring consistent and accurate metabolite quantification.

2.4. Comprehensive Transcriptome Sequencing and Data Analysis Workflow

RNA extraction was performed using the RNAprep Pure Plant Kit (Tiangen, Beijing, China) for plants and TRIzol Reagent (Life Technologies, Carisbad, CA, USA) for animals, followed by RNA quantification and integrity assessment using NanoDrop 2000 and the RNA Nano 6000 Assay Kit on the Agilent Bioanalyzer 2100 system. For transcriptome sequencing, 1 μg of RNA per sample was used to prepare sequencing libraries with the Hieff NGS Ultima Dual-mode mRNA Library Prep Kit, involving mRNA purification, cDNA synthesis, end-repair, adaptor ligation, PCR amplification, and quality assessment on the Agilent Bioanalyzer 2100 system. Libraries were sequenced on an Illumina NovaSeq platform to generate 150 bp paired-end reads. Raw reads were processed to obtain clean data through quality control steps, including adapter removal and low-quality read filtering, with mapping to the reference genome performed using Hisat2. Gene expression levels were quantified using FPKM, while differential expression analysis was conducted with DESeq2 (for biological replicates) or edgeR (for non-replicates) based on thresholds of adjusted p-value < 0.01 and Fold Change ≥ 2. KEGG pathway enrichment analysis was carried out using KOBAS and clusterProfiler, and alternative splicing events were quantified using rMATS software (version 4.3.0 turbo). This workflow ensured the generation of high-quality transcriptomic data and comprehensive analyses.

After sequencing, low-quality sequences were filtered out to ensure the quality of the aligned reads. Of the filtered sequences, 95% were successfully aligned, while 5% were unaligned.

2.5. Combined Analysis of Transcriptomic and Metabolomic Data

In the combined analysis of transcriptomic and metabolomic data, differentially expressed genes (DEGs) and differential metabolites (DMs) were identified using statistical thresholds (p < 0.05). DEGs and DMs were mapped to pathways through the KEGG database, with enzyme–gene associations established via EC numbers. Enriched pathways with significant overlaps of DEGs and DMs were identified using hypergeometric tests. Correlations between DEGs and DMs were calculated using Pearson correlation coefficients, retaining interaction pairs with r > 0.90 (p < 0.01). A network of these interactions was constructed and visualized using Cytoscape to reveal integrated molecular relationships.

The modules were defined as clusters of genes or metabolites that are co-expressed or correlated, as determined by Weighted Gene Co-expression Network Analysis (WGCNA). These modules represent coordinated patterns of biological processes and are key to understanding stress adaptation mechanisms.

2.6. qRT-PCR Validation

To validate the results of the differentially expressed genes (DEGs), 15 DEGs were randomly selected for measurement using qRT-PCR. Primers were designed using Primer Premier 5.0 software (

Table 1. qRT-PCR primer information.

NO.	Gene ID	Forward Primer	Reverse Primer
1	evm.TU.ctg3729.65	ATGGCTGCTTCTTCCAACAC	AGCTCATGAGACGTTCAATG
2	evm.TU.ctg8859.158	CAAAGACAGTTTGTGTCACG	TCTTGGATGGTGGCATGGAC
3	evm.TU.ctg10104.25	CTTCCTCAACTCCCTCACCT	GAAGTTGATTCGGCGTCGAT
4	evm.TU.ctg5993.28	AGAAGACAAGCTTTGGAGAG	CCCAATAGGGTACACTTTCT
5	evm.TU.ctg3507.25	TTTGGCCTCTGGTTTGGTTC	ACTGGTTCGGTGGCTACAAC
6	evm.TU.ctg6312.12	ACCCTCCAATTTCCAAAGTC	CTATGGTACCAATTACGAGG
7	evm.TU.ctg4592.109	ATGGGGGGTCTTTGTTCTAA	CCCTCCTTTACACTTGTCAA
8	evm.TU.ctg7070.226	CTTCACTTCACTCACTCACT	AGCGAGATGCTTTCAATGAG
9	evm.TU.ctg4437.300	ATGGGCACTGTGATTGACTC	CAAATTTGAGGAGTGCAGTG
10	evm.TU.ctg5993.15	CTGTTTTCAGGAGAGTTACC	GAACTCCCCGAAGTGTGAAA
11	evm.TU.ctg1707.227	CCATCTAACCAAACCCGACG	AACTCAACTTCGTCGTCAGG
12	evm.TU.ctg11098.136	GAGTTGTGCACTTAGATTGC	CGCTAGATGAGAATAGGAGA
13	evm.TU.ctg3686.26	TTGGTGGAGCTTTTTGTGAG	CTCCCTTTGCCAAATCCAAA
14	evm.TU.ctg10149.89	ATGGCGAGTAAAAGTGCTGA	TTTCCCGAGCTAGCATTGCC
15	evm.TU.ctg8474.31	GTTTTCTGAAAGATTGGGGG	CAACTTCTGGCTAACACTCC

). The qRT-PCR was performed following the instructions provided in the TB Green Premix Ex Taq™ II kit (Takara, Kusatsu, Japan).

The reaction system included 10 μL of TB Green Premix Ex Taq™ II, 1 μL of forward primer, 1 μL of reverse primer, 1 μL of cDNA, and 7 μL of RNase-Free H₂O, for a total volume of 20 μL. The relative expression levels of genes were calculated using the 2^−ΔΔCt method. All gene expression analyses were performed with three biological replicates [19].

3. Results

3.1. Metabolic Profile of Caucasian Clover at Three Time Points of Low-Temperature Stress

The analysis was conducted using the UPLC-ESI-MS/MS system (UPLC, Waters Acquity I-Class PLUS; MS, Applied Biosystems QTRAP 6500+), which identified 1150 metabolites across 12 samples (Supplementary Table S1). The metabolite classification analysis of Caucasian clover elucidated the key components of its metabolic network. According to the HMDB database classification, oxy-oxygenated compounds and carboxylic acids, along with their derivatives, were the predominant types of metabolites, comprising 71 and 70 metabolites, respectively. Furthermore, the quantities of fatty acyls, isoprenoid lipids, and flavonoid metabolites were notably high, with 58, 51, and 27 metabolites identified, respectively. In contrast, the number of phenolic compounds and phenol-like metabolites was relatively low (Figure 1a).

To further analyze the changes in metabolite composition, principal component analysis (PCA) was conducted to elucidate the overall trends in the metabolome under various low-temperature stress treatments. The results indicated that PC1 and PC2 accounted for 29.3% and 19.3% of the total variance, respectively. Groups A and B clustered in the upper left corner of the plot, reflecting similar early metabolic responses. In contrast, group C was positioned in the lower left corner, significantly separated from the control and early treatment groups, suggesting substantial metabolic changes during the mid-phase. Group D was distinctly located in the upper right corner, further emphasizing its unique metabolic regulation mechanism under low-temperature stress (Figure 1b).

Based on the metabolic differences identified through PCA, hierarchical clustering analysis of metabolite expression further elucidated the relationships among the samples. The clustering results demonstrated that samples A and B formed a tight cluster, indicating highly similar metabolic characteristics. In contrast, samples C and D emerged as independent branches, reflecting significant metabolic differentiation as the duration of low-temperature stress increased. The independent grouping of C and D also suggests that prolonged low-temperature stress induced more distinct and profound metabolic regulatory changes (Figure 1c).

Further K-means clustering analysis categorized the metabolites into two major clusters, comprising 548 and 772 metabolites, respectively. Metabolites in cluster 1 exhibited higher expression levels in groups A and B but demonstrated a significant decrease following treatments C and D under low-temperature stress. In contrast, metabolites in cluster 2 displayed more stable overall expression across all treatment groups, with only minor fluctuations at certain time points and no significant trend changes (Figure 1d).

3.2. Identification of Differentially Accumulated Metabolites at Different Stages of Low-Temperature Stress

To identify key differentially accumulated metabolites at various stages, supervised orthogonal partial least-squares discriminant analysis (OPLS-DA) was employed to construct a classification model and differentiate metabolites across these stages. The OPLS-DA revealed significant changes in the metabolic profiles of Caucasian clover under different stages of low-temperature stress. Notable separation was observed between the control group and the treatment groups, as well as among the treatment groups themselves. Furthermore, as the duration of low-temperature stress increased, the metabolic changes became increasingly pronounced. The explanation rate of the first principal component progressively increased from 42% to 49%, indicating that metabolic changes accumulated over time. Biological replicates within each group were tightly clustered, ensuring the reliability of the results (Figure 2a–c).

To quantify the dynamic changes in metabolite expression, we systematically compared the differences in metabolite expression across various durations of low-temperature stress. As the treatment duration increased, the number of upregulated metabolites rose from 150 to 217, while the number of downregulated metabolites also increased from 223 to 244 (Figure 2d).

3.3. KEGG Pathway Enrichment Analysis of Differential Metabolites

To investigate the dynamic variations in differential metabolites under various low-temperature stress conditions, we analyzed the distribution of upregulated and downregulated metabolites, as well as their enrichment in metabolic pathways. This analysis aims to elucidate the stage-specific response characteristics of Caucasian clover to low-temperature stress.

Venn diagrams illustrating the distribution and characteristics of upregulated and downregulated metabolites under different low-temperature stress conditions were utilized. In terms of upregulated metabolites, we identified 100 specific metabolites in the A vs. B comparison, 119 in B vs. C, and the highest number, 165, in C vs. D. Furthermore, overlaps between comparison groups were observed. For instance, 25 metabolites consistently exhibited upregulation in both the A vs. B and B vs. C comparisons, with only one metabolite showing consistent upregulation across all three comparisons (Figure 3a). The analysis of downregulated metabolites revealed 139 specific metabolites in the A vs. B comparison, 125 in B vs. C, and 155 in C vs. D. The overlaps among comparison groups further indicated shared responses of specific metabolites. Notably, 36 metabolites were downregulated in both the A vs. B and B vs. C comparisons, while 8 metabolites demonstrated consistent downregulation across all three comparisons (Figure 3b).

Further KEGG pathway enrichment analysis provided deeper insights into the functional roles of the differential metabolites. In the comparison of A vs. B, differential metabolites were significantly enriched in C5-branched dibasic acid metabolism, alpha-linolenic acid metabolism, arginine and proline metabolism, as well as glyoxylic acid and dicarboxylic acid metabolism. As the treatment duration increased, the comparison of B vs. C revealed significant enrichment in isoflavonoid biosynthesis, phenylalanine metabolism, and purine metabolism. In the comparison of C vs. D, metabolites were significantly enriched in C5-branched dibasic acid metabolism, alpha-linolenic acid metabolism, and glyoxylic acid and dicarboxylic acid metabolism (Figure 3c–e).

3.4. Transcriptomics Analysis of Caucasian Clover Under Different Low-Temperature Treatments

3.4.1. Analysis of Differentially Expressed Genes

The transcriptomic response of Caucasian clover to various low-temperature stress conditions demonstrated significant alterations in gene expression. In the comparison of A vs. B, 12,401 genes were found to be upregulated, while 12,234 genes were downregulated. As the duration of low-temperature treatment increased, in the comparison of B vs. C, the number of upregulated genes decreased to 6105, whereas the number of downregulated genes increased to 6506. In the comparison of C vs. D, the number of upregulated genes further declined to 4991, and the number of downregulated genes rose to 6511 (Figure 4a).

Among the upregulated genes, distinct shared and specific characteristics were observed across the comparison groups. In the comparison of A vs. B, 10,999 specific genes were upregulated; in B vs. C, 5462; and in C vs. D, 4197. Additionally, 537 genes were upregulated in both A vs. B and B vs. C, while 10 genes consistently exhibited upregulation across all three comparisons (Figure 4b). Regarding downregulated genes, 10,852 specific genes were downregulated in A vs. B, 5541 in B vs. C, and 4895 in C vs. D. Among the shared downregulated genes, 60 showed consistent downregulation across all three comparisons, and 419 genes were downregulated in both A vs. B and B vs. C (Figure 4c).

3.4.2. KEGG Pathway Enrichment Analysis of Differentially Expressed Genes

The KEGG enrichment analysis of differentially expressed genes revealed dynamic changes in the metabolic and biological pathways of Caucasian clover during various stages of low-temperature stress. In the comparison of A vs. B, the significantly enriched pathways included circadian rhythm–plant, photosynthesis–antenna proteins, alpha-linolenic acid metabolism, and monoterpenoid biosynthesis (Figure 5a).

As the duration of low-temperature stress extended to 6 h (B vs. C), the alterations in plant metabolic and biological pathways became more pronounced. At this stage, the significantly enriched pathways included circadian rhythm–plant, alpha-linolenic acid metabolism, photosynthesis–antenna proteins, and linoleic acid metabolism (Figure 5b). When the duration of low-temperature stress was extended to 12 h (C vs. D), further adjustments in metabolic and regulatory pathways were observed, with significantly enriched pathways including photosynthesis–antenna proteins, circadian rhythm–plant, photosynthesis, and linoleic acid metabolism (Figure 5c).

3.4.3. Validation by Quantitative Real-Time PCR (qRT-PCR)

This study randomly selected 15 differentially expressed genes to validate the impact of cold stress on key gene expression in Caucasian clover using qRT-PCR. The results are shown in Figure 6. The expression dynamics of these genes at various stages of cold stress suggest their potential involvement in cold adaptation. For example, the expression of genes evm.TU.ctg3729.65 and evm.TU.ctg10149.89 increased rapidly at the early stress stage (2 h), potentially playing a key role in signal transduction. In contrast, genes evm.TU.ctg4592.109 and evm.TU.ctg8474.31 peaked in expression at the later stage of stress (12 h), indicating their role in long-term stress adaptation. These validated genes are involved in several biological processes, including antioxidant defense, membrane lipid repair, and signal transduction, and they align with KEGG pathway analysis results. For instance, the expression of evm.TU.ctg6312.12 is strongly associated with the activity of the ABA signaling pathway and may enhance plant cold resistance by regulating osmotic adjustment and antioxidant protection.

Furthermore, the dynamic changes in these validated genes further support the conclusions of the integrated transcriptomic and metabolomic analysis, indicating that the accumulation of antioxidant metabolites and the coordinated expression of related genes are central to the cold stress resistance mechanisms in Caucasian clover.

3.5. Integrated KEGG Pathway Enrichment Analysis of Transcriptomics and Metabolomics

The integrated KEGG pathway enrichment analysis of transcriptomics and metabolomics at various stages of low-temperature stress in Caucasian clover reveals a multi-layered metabolic regulatory mechanism, transitioning from a rapid response to gradual adaptation. During the early stage (A vs. B), the plant primarily activates alpha-linolenic acid metabolism along with arginine and proline metabolism, thereby adjusting membrane lipid composition and accumulating osmoregulatory compounds to maintain cellular stability. Furthermore, secondary metabolic pathways, including phenylpropanoid and flavonoid biosynthesis, are activated to scavenge reactive oxygen species (ROS) and protect cells from oxidative damage. The enrichment of carbon metabolism and carbon fixation pathways in photosynthetic organisms indicates that the plant sustains growth and fundamental metabolic activities by regulating energy metabolism (Figure 7a).

During the mid-phase (B vs. C), the metabolic network exhibits more intricate regulation. Lipid metabolism, including linoleic acid and alpha-linolenic acid metabolism, remains enriched, thereby enhancing membrane fluidity and stability. Additionally, glutathione metabolism is significantly upregulated, indicating the sustained activation of the plant’s antioxidant defense system. The enrichment of carbohydrate metabolism, encompassing starch and sucrose metabolism as well as amino sugar and nucleotide sugar metabolism, along with amino acid metabolism, such as valine, leucine, and isoleucine degradation and arginine and proline metabolism, reflects the plant’s strategy to manage stress by supplying energy and regulating cellular osmotic pressure. In the mid-phase, the activation of terpenoid metabolism, including terpenoid backbone biosynthesis and sesquiterpenoid and triterpenoid biosynthesis, alongside signal transduction pathways, particularly plant hormone signal transduction involving ABA and JA signaling, was observed. This suggests that the plant modulates the expression of stress-responsive genes through signaling molecules (Figure 7b).

In the later phase (C vs. D), the plant enters a stable adaptation phase characterized by the ongoing enrichment of energy metabolic pathways, including the tricarboxylic acid (TCA) cycle and carbon fixation processes, which are vital for sustaining energy supply during extended periods of stress. Simultaneously, secondary metabolic pathways, such as phenylpropanoid and flavonoid biosynthesis, are further enhanced to support reactive oxygen species (ROS) detoxification and to strengthen the cell wall. The persistent activity of glutathione metabolism highlights the essential role of the antioxidant system in adapting to long-term low temperatures. Furthermore, the continuous enhancement of arginine and proline metabolism underscores the multifaceted roles of proline in osmotic protection, antioxidation, and energy storage (Figure 7c).

3.5.1. Weighted Gene Co-Expression Network Analysis (WGCNA) of Key Metabolites and Gene Co-Expression Networks

The gene clustering dendrogram generated through Weighted Gene Co-expression Network Analysis (WGCNA) clarified gene co-expression patterns by categorizing them into multiple modules, each displaying distinct expression profiles under low-temperature stress conditions. Dynamic tree cutting generated an initial set of modules, reflecting the primary clustering relationships among genes. Subsequently, modules were merged based on the similarity of their eigengenes, resulting in a set of biologically meaningful modules. Color-distinguished modules indicate variations in gene expression patterns. Specifically, the blue module shows a strong correlation with low-temperature stress, indicating its pivotal role in antioxidant defense and signal transduction pathways. The green module is significantly associated with metabolites such as jasmonic acid and L-proline, suggesting its crucial involvement in lipid metabolism and osmotic regulation. Additionally, the red module is involved in hormone regulatory pathways, including abscisic acid (ABA) and jasmonic acid (JA) (Figure 8a). Further analysis revealed that the eigengenes of distinct modules were significantly correlated with low-temperature stress and metabolite accumulation, suggesting that the genes within these modules collaboratively orchestrate the plant’s adaptive mechanisms to low-temperature stress.

By utilizing gene clustering dendrograms and heatmap analyses, the integrated expression patterns of selected metabolites and genes elucidated the co-expression of modules associated with metabolites and their distinct expression profiles under low-temperature stress conditions. The clustering dendrogram demonstrated that genes were divided into multiple modules based on expression similarity. Color-coded modules reflect the potential functional associations of distinct gene sets in metabolite regulation. Heatmap analysis revealed that gene modules associated with metabolites, including L-ascorbic acid, jasmonic acid, and L-proline, exhibited significant differences in expression patterns under stress conditions, with red indicating high expression and blue indicating low expression. Notably, modules significantly associated with jasmonic acid and L-proline, specifically the green and blue modules, exhibited pronounced stress responses, suggesting their potential roles in antioxidant defense, lipid metabolism, and signal transduction pathways. Furthermore, modules with high expression levels related to L-ascorbic acid may be associated with the antioxidant function of reactive oxygen species (ROS) scavenging (Figure 8b). Comprehensive analysis indicates that genes within these metabolite-associated modules play critical roles in the molecular adaptation mechanisms of plants under low-temperature stress.

3.5.2. Heatmap Analysis of Correlation Between Module Eigengenes and Metabolites

Using WGCNA, we identified correlations between multiple gene modules and specific metabolites, elucidating the molecular adaptation mechanisms of Caucasian clover under low-temperature stress conditions. The MEorangered4 module showed a correlation coefficient of 0.64 with ascorbic acid (AsA), indicating its potential role in antioxidant defense mechanisms that help plants mitigate oxidative stress induced by low-temperature conditions. The MEplum1 and MEbrown modules demonstrated correlation coefficients of 0.78 and 0.64 with D-proline (D-Pro), respectively, highlighting proline’s pivotal role in osmotic regulation and antioxidant defense mechanisms. The MEyellowgreen and MEturquoise modules showed correlation coefficients of 0.67 and 0.61 with abscisic acid (ABA), respectively, indicating their involvement in ABA signaling pathways and the regulation of water balance and osmotic adjustments in plants under low-temperature stress. The MEdarkorange module exhibited a high correlation coefficient of 0.86 with Maltose, suggesting that this module may facilitate plant adaptation to low-temperature stress through energy metabolism and osmotic regulation mechanisms. The MEgrey60 module showed a correlation coefficient of 0.91 with Diphospho-D-Glucose, highlighting its critical role in energy storage and metabolic regulation, especially in maintaining cellular energy supply under low-temperature stress. The MElightsteelblue1 module showed a correlation coefficient of 0.97 with D-Fructose, indicating that fructose plays a central role in sugar metabolism and osmotic regulation under low-temperature stress. The MEroyalblue module showed a correlation coefficient of 0.68 with Maltose, further substantiating the role of carbohydrate metabolism in adaptation to low-temperature stress. The MEdarkolivegreen module showed a correlation coefficient of 0.63 with jasmonic acid (JA), indicating its potential involvement in the low-temperature stress response through signal transduction pathways. The MEblack module showed a correlation coefficient of 0.78 with Glucose 6-Phosphate (G6P), revealing its significant role in sugar metabolism and energy regulation. Finally, the MEpink module showed a high correlation coefficient of 0.94 with 3-O-Galloyl-D-glucose, suggesting that this module may enhance plant cold resistance through antioxidant activity and the synthesis of secondary metabolites (Figure 9). Based on the aforementioned analysis, this study selected the MEpink and MEplum1 modules for subsequent comprehensive analysis.

3.5.3. KEGG Enrichment Analysis of Key Modules and Construction of Core Gene Interaction Networks

Using KEGG enrichment analysis of genes within the Pink module, we identified multiple metabolic pathways essential for plant responses to low-temperature stress. The enrichment of these pathways elucidates how plants enhance their adaptability to low temperatures by modulating molecular and metabolic network activities. First, the significant enrichment of the glycerophospholipid metabolism and sphingolipid metabolism pathways suggests that low-temperature stress may enhance cell membrane stability by modulating the composition and fluidity of membrane lipids, thereby protecting cells from cold-induced damage. The modulation of membrane lipid metabolism helps maintain the integrity of cellular structures and ensures the functional stability of cell membranes under low-temperature conditions. Concurrently, the enrichment of the starch and sucrose metabolism pathways, as well as the amino sugar and nucleotide sugar metabolism pathways, indicates that plants regulate carbon metabolism and carbohydrate synthesis to provide essential energy support, thereby mitigating the metabolic stress induced by low-temperature conditions. The enrichment of the glutathione metabolism pathway further underscores the antioxidant mechanisms of plants, suggesting that plants enhance antioxidant reactions to eliminate reactive oxygen species (ROS) induced by low temperatures, thereby protecting cells from oxidative damage (Figure 10a). Further gene interaction network analysis identified two pivotal genes, EVM-20prediction-20ctg2798.293 and EVM-20prediction-20ctg1205.46, which showed high connectivity within the network. This high connectivity indicates their central roles in the response to low-temperature stress, suggesting they may serve as crucial regulatory factors in plant adaptation to cold conditions. EVM-20prediction-20ctg2798.293 may enhance plant adaptation to low temperatures by modulating protein synthesis and signal transduction pathways. The high connectivity within the network indicates interactions with multiple other metabolic pathways, potentially facilitating protein synthesis and translation under low-temperature conditions, thereby aiding plants in maintaining cellular functional stability. Conversely, EVM-20prediction-20ctg1205.46, a multidrug resistance-associated protein of the MATE family, may aid plants in eliminating harmful substances and reducing oxidative damage induced by low temperatures by participating in detoxification and antioxidant reactions under cold conditions. The synergistic interaction of the two genes may enhance plant tolerance to low temperatures by modulating protein synthesis, antioxidant responses, and metabolic pathways. The expression patterns of these core genes under low-temperature stress were further validated through heatmap analysis (Figure 10c). Heatmap analysis revealed that the expression levels of both genes gradually increased with prolonged low-temperature treatment, reaching their peak at 12 h of low-temperature exposure (Group D). This suggests that these two genes may play significant roles in long-term adaptation to low temperature. Particularly, in Group B (2 h) and Group C (6 h) of low-temperature exposure, the expression levels of EVM-20prediction-20ctg2798.293 gradually increased, suggesting its involvement in the rapid response of plants to low temperature during the initial and mid-term phases. Conversely, EVM-20prediction-20ctg1205.46 exhibited a progressively increasing trend across all treatment groups, peaking at 12 h of low-temperature exposure (Group D). This indicates that it may aid plants in coping with cold damage during long-term adaptation to low-temperature stress by modulating cellular repair and antioxidant responses (Figure 10e).

KEGG enrichment analysis of the Plum module revealed multiple key metabolic pathways, particularly those involved in lipid metabolism, antioxidant responses, and signal transduction. Specifically, the significant enrichment of alpha-linolenic acid metabolism and linoleic acid metabolism pathways suggests that the Plum module may enhance the fluidity and stability of cell membranes by modulating fatty acid metabolism, thereby mitigating damage caused by low-temperature stress. Additionally, the enrichment of the MAPK signaling pathway underscores the role of module genes in signal transduction during low-temperature stress responses, indicating that plants may regulate the expression of downstream genes through the MAPK signaling pathway to initiate stress resistance mechanisms. The enrichment of the glutathione metabolism and peroxisome-related pathways indicates that the module plays a crucial role in scavenging reactive oxygen species (ROS) and alleviating oxidative damage induced by low temperatures. The enrichment results of these pathways indicate that the Plum module enhances plant adaptability to low-temperature stress through the synergistic interaction of multiple metabolic and signaling pathways (Figure 10b). Further gene interaction network analysis identified the core gene within the Plum module: EVM-20prediction-20ctg4298.8. Its high connectivity within the network suggests that this gene may play a pivotal regulatory role in low-temperature responses. Specifically, EVM-20prediction-20ctg1327.25, a gene involved in the MAPK signaling pathway, and EVM-20prediction-20ctg87.24, associated with fatty acid metabolism, may maintain cellular stability under low-temperature stress by modulating membrane lipid composition and signal transduction (Figure 10d). In the gene expression heatmap analysis of the Plum module, expression changes of genes at different low-temperature treatment time points were observed. In Group B (2 h of low-temperature exposure), the expression levels of multiple pivotal genes significantly increased, indicating that these genes may participate in the initial response mechanisms to low temperature. These genes may aid plants in coping with short-term low-temperature stress by rapidly initiating signal transduction, antioxidant responses, and protein synthesis. In Group C (6 h) and Group D (12 h) of low-temperature exposure, although the expression levels of some genes remained high, there was no further significant increase. This indicates that these genes may have already initiated adaptation mechanisms during the initial phase of low temperature, with the later phase primarily maintaining these response processes (Figure 10f). This dynamic change reveals the time-dependent nature of gene expression under low-temperature stress, wherein certain genes play pivotal regulatory roles during the early stages of stress.

3.6. DAMs and DEGs Involved in Isoflavonoid Biosynthesis Pathway

Based on the integrated analysis of the transcriptome and metabolome, we concluded that flavonoid metabolism, particularly isoflavonoid biosynthesis, plays a critical role in the low-temperature stress tolerance of Caucasian clover. Accordingly, we delineated an isoflavonoid biosynthesis pathway that incorporates key metabolites and genes, utilizing the KEGG database, which includes 11 metabolites and 461 genes.

Under low-temperature stress, the isoflavonoid metabolism pathway exhibits a highly dynamic regulatory pattern, with gene expression and metabolite accumulation undergoing stage-specific changes that reflect the molecular mechanisms involved in the plant’s response to low-temperature stress. During the initial stage of stress (2 h), the plant rapidly activates key upstream genes in the pathway, particularly Isoflavone Synthase (IFS), whose significant upregulation triggers the core processes of isoflavonoid metabolism. This activation leads to the generation of substantial amounts of precursors such as Daidzein and Genistein, thereby facilitating the accumulation of downstream metabolites. Concurrently, the expression of Isoflavone-4′-O-Methyltransferase (HI4OMT) and Isoflavone-7-O-Glucosyltransferase (IF7GT) is markedly upregulated, promoting the production of more stable and biologically active isoflavonoid derivatives, including Formononetin (methoxy-Daidzein) and Formononetin 7-O-glucoside, through methylation and glycosylation modifications.

The activation of downstream enzymes, such as CYP81E1 and VR (Vestitone Reductase), enhances the synthesis of the antimicrobial metabolite Medicarpin, thereby strengthening the plant’s defense mechanisms. At this stage, metabolites like Glycitin and Formononetin 7-O-glucoside reach their peak concentrations, indicating that the plant improves the stability and functional durability of these metabolites through glycosylation and modification processes. Concurrently, the accumulation of Medicarpin increases significantly, suggesting its role in enhancing the plant’s adaptive capacity during the mid-phase of stress through antimicrobial and protective functions. Under prolonged stress (12 h), the pathway exhibits a different regulatory pattern: the expression of certain genes, such as CYP81E1, peaks, indicating the plant’s prioritization of antimicrobial metabolite accumulation to counter pathogen invasion. In contrast, the expression of genes like HIDH and VR significantly declines, possibly reflecting the plant’s reallocation of resources toward other stress-resistance pathways. During this period, the concentration of Medicarpin continues to rise, serving as the core protective metabolite for the plant’s long-term stress resistance, while the sustained accumulation of metabolites such as Biochanin A and Formononetin likely contributes to maintaining cellular homeostasis through antioxidant and signaling regulation. However, the concentration of Glycitin declines, indicating a reduced reliance on it by the plant (Figure 11).

4. Discussion

In the context of global climate change, low-temperature stress has emerged as a significant abiotic stressor affecting plant growth, development, and yield [20,21]. To adapt to low-temperature stress, plants utilize complex metabolic regulatory networks that involve the accumulation of protective metabolites, modulation of antioxidant defenses, and activation of secondary metabolic pathways [22,23]. Caucasian clover, a crucial forage resource, demonstrates metabolic response characteristics under low-temperature stress that are vital for understanding its adaptability. This knowledge also provides a theoretical foundation for exploring the molecular mechanisms of low-temperature tolerance in plants and informs related breeding strategies. However, research on the dynamic metabolic response of Caucasian clover to low-temperature stress remains relatively limited, highlighting a valuable opportunity for the current investigation.

4.1. The Response of DAMs in Caucasian Clover to Low-Temperature Stress

This study reveals significant metabolic dynamics and regulatory features of Caucasian clover under low-temperature stress, demonstrating a stage-specific metabolic network response pattern. By integrating metabolite classification, principal component analysis (PCA), K-means clustering, and KEGG enrichment analysis, the metabolic strategies and adaptive mechanisms of the plant in response to low-temperature stress are delineated, providing new insights into the molecular foundation of plant responses to such stress. The metabolite classification analysis indicates that oxygenated compounds and carboxylic acids, along with their derivatives, are predominant, which correlates with increased energy metabolism and antioxidant demand. This suggests that the plant prioritizes metabolic protective mechanisms under low-temperature stress. This observation aligns with the critical role of carboxylic acid metabolism in low-temperature adaptation, as highlighted in studies of Arachis hypogaea [24] and Brassica napus [25], emphasizing the conserved nature of core metabolite classification functions across different plant species. Additionally, the elevated levels of fatty acyl metabolites further suggest that Caucasian clover effectively mitigates cold-induced damage to cellular structures by modulating membrane lipid composition, thereby preserving membrane stability. This finding is consistent with the membrane protection mechanisms observed in studies of Actinidia arguta [21] and Nicotiana tabacum [10].

PCA analysis and heatmap clustering results elucidated the metabolic dynamics of Caucasian clover across various stages of low-temperature stress. During the short-term stress phase (A vs. B), the metabolite composition remained relatively consistent, likely reflecting an initial adaptation characterized by the rapid accumulation of soluble sugars and amino acids to maintain osmotic pressure and support fundamental metabolic functions. This swift response aligns with findings in wheat, where adaptation to low-temperature stress was achieved promptly through proline and sugar metabolism [26]. As the duration of stress increased, the metabolite distribution in the mid-phase (C) and late-phase (D) exhibited significant differentiation, indicating a transition from a rapid response to more complex metabolic network regulation. In this latter phase, secondary metabolic pathways, including those related to phenylpropanoids and flavonoids, were notably activated, providing essential metabolic support for reactive oxygen species (ROS) scavenging and enhanced antioxidant defense [22]. This trend was corroborated by studies on Arachis hypogaea, where the activation of secondary metabolic pathways became particularly evident after 6 h [24]. K-means clustering analysis further clarified the functional roles of metabolites at various stages. Metabolites in Cluster 1 were highly expressed during the early stress phase but significantly decreased in later stages, suggesting their involvement in short-term adaptive functions such as energy provision and osmotic regulation [27]. In contrast, metabolites in Cluster 2 were consistently expressed throughout the stress period, likely playing a crucial role in sustaining long-term adaptation and basal metabolic functions. This functional differentiation is consistent with findings in Medicago sativa [28], Brassica napus [29], and Arachis hypogaea [24], indicating that the stable accumulation of secondary metabolites is vital for stress tolerance [18].

The dynamic changes in metabolite quantities further substantiate this regulatory pattern. As the duration of low-temperature treatment increased, the number of upregulated metabolites rose from 150 to 217, while the number of downregulated metabolites increased from 223 to 244. This trend indicates that Caucasian clover activates adaptive regulatory pathways by upregulating metabolites in response to low-temperature stress [30], while simultaneously downregulating specific metabolites to optimize energy allocation and prioritize essential metabolic processes [10,31]. This phenomenon is further supported by research on wheat, where the sustained accumulation of proline and soluble sugars is recognized as a key mechanism for mitigating oxidative stress [26]. The significant separation observed in OPLS-DA analysis further elucidated the time-dependent metabolic characteristics at different stress stages. In the early stages, rapid regulation of metabolites predominated, enabling the plant to swiftly adapt to environmental changes. In contrast, during the middle and late stages, the accumulation of secondary metabolites enhanced antioxidant defense and signal transduction, providing critical support for long-term stress adaptation. This stage-dependent regulatory pattern is not only evident in Caucasian clover but has also been widely validated in cold-tolerant plants such as Vicia sativa [30], Brassica napus [32], and Medicago sativa [28], suggesting that it may represent a universal low-temperature adaptation strategy among plants.

KEGG pathway enrichment analysis further elucidated the time-dependent regulatory patterns of differential metabolites. In the early stages of low-temperature stress (A vs. B), differential metabolites were significantly enriched in pathways such as C5-branched dicarboxylic acid metabolism, alpha-linolenic acid metabolism, arginine and proline metabolism, and glycolate and dicarboxylic acid metabolism. The enrichment of these pathways primarily reflects the plant’s rapid adaptive response to acute low-temperature stress. The significant enrichment of arginine and proline metabolism suggests that proline may serve as a key osmoregulatory compound under low-temperature stress, stabilizing proteins and membrane structures, scavenging reactive oxygen species (ROS), and regulating intracellular osmotic pressure. This phenomenon has been extensively documented in plants such as Actinidia arguta [21] and Triticum aestivum [26], demonstrating that proline accumulation is a conserved strategy for plants to rapidly adapt to low-temperature stress. Additionally, the activity of alpha-linolenic acid metabolism underscores the critical role of membrane lipid remodeling under low temperatures, as the plant maintains membrane stability by adjusting lipid composition, thereby mitigating damage to cellular structures induced by low-temperature stress. This finding is consistent with studies on Cocos nucifera [33] and Nicotiana tabacum [10] under low-temperature stress, indicating that the activation of membrane lipid metabolism is a key adaptive response of plants to low-temperature stress. In the middle stage (B vs. C), the enrichment of differential metabolites revealed a significant functional shift, with isoflavonoid biosynthesis, phenylalanine metabolism, and purine metabolism identified as the major enriched pathways. The pronounced enhancement of isoflavonoid biosynthesis and phenylalanine metabolism indicates that the plant increases its antioxidant capacity through the accumulation of secondary metabolites, such as flavonoids and phenylpropanoids. These metabolites not only scavenge reactive oxygen species (ROS) but also strengthen the cell wall structure, thereby improving the plant’s tolerance to low-temperature stress. This mechanism has been extensively validated in studies on Arachis hypogaea [24] and Brassica napus [34], further emphasizing the crucial role of secondary metabolic pathways in cold resistance. Moreover, the significant enrichment of purine metabolism may be associated with the heightened energy metabolism demands of the plant during the mid-stress phase. Purine metabolism is intricately linked to intracellular nucleic acid and energy metabolism, providing both material and energy support for stress adaptation, particularly in cell repair and signal transduction processes. In the later stages of low-temperature stress (C vs. D), differential metabolites were significantly enriched in the C5-branched dicarboxylic acid metabolism, alpha-linolenic acid metabolism, and glycolate and dicarboxylic acid metabolism pathways. Unlike the early stage, the sustained activity of these pathways in the later stages underscores the importance of core metabolic processes in maintaining cellular homeostasis and energy supply. The significant enrichment of these core metabolic pathways indicates that plants enhance their adaptation to prolonged low-temperature stress by optimizing resource allocation to support essential physiological processes, such as the sustained expression of antioxidant enzymes and cellular repair. Furthermore, the sustained activity of glycolate and dicarboxylic acid metabolism may be closely linked to the reallocation of carbon metabolism and the maintenance of oxidative balance. This observation aligns with findings reported in studies on Phaseolus vulgaris [35] and Cucumis sativus [36] under low-temperature stress, further confirming the conserved role of core metabolic pathways in long-term cold adaptation.

4.2. The Response of DEGs in Caucasian Clover to Low-Temperature Stress

The transcriptomic response of Caucasian clover to low-temperature stress revealed significant changes in gene expression. As the duration of low-temperature treatment increased, the patterns of upregulation and downregulation exhibited distinct phase-dependent variations. In the comparison of A vs. B, 12,401 genes were upregulated, while 12,234 genes were downregulated, indicating a rapid response characterized by the activation of numerous genes to counteract low-temperature stress. With the extension of low-temperature treatment, in the comparison of B vs. C, the number of upregulated genes decreased to 6105, whereas the number of downregulated genes increased to 6506, reflecting a gradual stabilization of the plant’s adaptation and an entry into a phase of relative balance. During the 12 h low-temperature treatment phase (C vs. D), the number of upregulated genes further decreased to 4991, while downregulated genes increased to 6511, suggesting that the plant’s response to prolonged low-temperature stress shifted towards resource redistribution and adaptive regulation.

Analysis of the common and specific upregulated genes across different comparison groups revealed that, during the early stage, 10,999 specific upregulated genes were significantly activated, with this number gradually decreasing as low-temperature treatment progressed. This observation suggests that, in the initial phase of low-temperature stress, plants must activate a substantial number of defense genes, particularly those involved in antioxidation, photosynthesis, membrane lipid remodeling, and osmoregulation [37]. This finding aligns with studies on crops such as Vigna radiata [38] and Cicer arietinum [39], indicating that during the early stages of low-temperature stress, plants maintain cellular stability by enhancing the antioxidative system, increasing osmotic protectants, and regulating membrane lipid metabolism. Under prolonged low-temperature stress, plants progressively optimize these responses, thereby reducing energy expenditure [40]. As plants transition into the adaptation phase, the activation of defense genes gradually declines, while key genes responsible for regulating basal metabolism and sustaining growth are activated to stabilize the physiological state [41]. Furthermore, 537 genes were persistently upregulated in comparisons A vs. B and B vs. C, suggesting that these genes likely play a role in the cold tolerance response during both the early and mid-phases of low-temperature stress, demonstrating robust cold resistance and adaptability. The downregulation pattern also exhibited a time-dependent trend. In the comparison of A vs. B, a total of 10,852 specific downregulated genes were identified, indicating the plant’s adaptive strategy to low-temperature stress by suppressing non-essential metabolic pathways and activating energy-saving mechanisms during the early phase [42]. As the low-temperature treatment progressed (B vs. C and C vs. D), certain genes remained downregulated, particularly those associated with growth and development, including genes involved in jasmonic acid and abscisic acid signaling pathways. This observation suggests that plants reduce energy consumption under prolonged low-temperature stress, prioritizing the maintenance of essential functions [43].

KEGG pathway enrichment analysis of differentially expressed genes revealed dynamic shifts in metabolic and biological pathways in Caucasian clover under low-temperature stress at various stages. In the early phase of low-temperature stress, the plant initially adjusts its physiological rhythm by activating genes associated with the circadian clock, thereby facilitating rapid adaptation to environmental fluctuations [44]. This mechanism serves as a central strategy for plants to maintain internal homeostasis in response to external environmental changes. Concurrently, the enrichment of photosynthesis–antenna proteins indicates that the plant optimizes energy capture and conversion by regulating the photosynthetic system during the early phase of low-temperature stress, ensuring an adequate energy supply under these conditions [45]. Although low-temperature stress generally suppresses photosynthesis, the plant finely tunes photosynthetic pathways in the initial phase to sustain essential metabolic processes.

The activation of lipid metabolism pathways, particularly the metabolism of alpha-linolenic acid, plays a critical role under low-temperature stress. Research indicates that low-temperature stress modulates lipid metabolism, leading to alterations in membrane lipid composition that enhance both membrane fluidity and stability, thereby mitigating cold-induced membrane damage [46,47]. This mechanism enables plants to maintain membrane integrity under stress conditions, thus preventing freeze–thaw damage. Furthermore, the enrichment of monoterpenoid biosynthesis suggests that plants activate secondary metabolic pathways in the early phases of low-temperature stress, producing compounds with antioxidant and antimicrobial properties that are essential for bolstering the plant’s defense against oxidative stress and pathogen attacks induced by low temperatures [48].

As the duration of low-temperature treatment increases, particularly from 2 to 6 h, the activation of metabolic pathways undergoes significant adjustments. While the plant circadian rhythm and photosynthesis–antenna protein pathways remain enriched, linoleic acid metabolism emerges as a crucial regulatory pathway over time. Linoleic acid is a vital polyunsaturated fatty acid in plants, playing a key role in maintaining membrane stability and modulating the antioxidant system [49]. This shift indicates the plant’s ongoing optimization of energy and lipid metabolism under prolonged stress, which contributes to enhanced low-temperature tolerance and stress resilience.

During the prolonged low-temperature stress phase (12 h), photosynthesis-related genes remain highly expressed. However, the enrichment of metabolic pathways indicates that the plant gradually transitions into a long-term adaptation mode, particularly characterized by the sustained enrichment of photosynthesis and linoleic acid metabolism. This suggests that the plant stabilizes cellular function by maintaining photosynthetic efficiency and optimizing lipid metabolism throughout this phase [50]. This regulatory mechanism closely resembles the low-temperature adaptation strategies observed in Vicia sativa [51] and Anthurium andraeanum [52], illustrating how plants enhance physiological adaptability and survival through coordinated metabolic reprogramming in response to low-temperature stress.

4.3. WGCNA Reveals the Low-Temperature Adaptation Mechanisms of Caucasian Clover

This study elucidates the molecular responses of Caucasian clover to low-temperature stress by integrating wide-targeted metabolomics and transcriptomics, using Weighted Gene Co-expression Network Analysis (WGCNA). WGCNA categorized all genes into multiple co-expression modules, with each significantly associated with antioxidant defense, lipid metabolism, and hormone signal transduction. This demonstrates the pivotal roles of these biological processes in plant adaptation to low temperatures. For example, the strong correlation between the MEorangered4 module and antioxidant pathways supports the idea that plants protect cellular structures under low-temperature stress by enhancing their antioxidant capacity to eliminate reactive oxygen species (ROS). This finding is consistent with studies underscoring the importance of antioxidant enzyme systems in plant cold resistance [53]. The strong correlation of the MEplum1 module with proline and jasmonic acid (JA) highlights proline’s dual role in osmotic regulation and antioxidant defense. This aligns with proline’s known role as an osmoprotectant in other plant species [50] and also corroborates JA’s function as a stress-related plant hormone. Furthermore, the correlation between the MEyellowgreen module and abscisic acid (ABA) indicates that ABA signal transduction critically regulates stomatal closure and osmotic adjustment in plants. This finding aligns with ABA’s central regulatory role in plant stress resistance [54].

The high correlation (0.94) between the Pink module and 3-O-Galloyl-D-glucose suggests that this module is pivotal to low-temperature adaptation in Caucasian clover. 3-O-Galloyl-D-glucose is a secondary metabolite with strong antioxidant activity. Its accumulation eliminates low-temperature-induced ROS, thereby protecting cells from oxidative damage. This finding aligns with the known role of secondary metabolites in plant stress resistance. For instance, catechin compounds in Arabidopsis significantly contribute to antioxidant defense mechanisms [55]. Furthermore, KEGG enrichment of the Pink module revealed significant enrichment in the glycerophospholipid, sphingolipid, starch and sucrose, and glutathione metabolism pathways. The active involvement of these pathways shows how plants enhance cell membrane fluidity and stability by adjusting lipid composition and energy metabolism. This prevents membrane rupture and functional disorders caused by low temperatures [50]. Regulating starch and sucrose metabolism provides essential energy and maintains cellular osmotic balance, thereby enhancing plant survival in low-temperature environments [51]. Glutathione metabolism enrichment reinforces antioxidant mechanisms by removing ROS, reducing oxidative damage, and enhancing overall stress resistance [39]. Gene interaction network analysis identified the core genes EVM-20prediction-20ctg2798.293 and EVM-20prediction-20ctg1205.46 within the Pink module. The former may enhance cold resistance by modulating protein synthesis and signal transduction, while the latter, a MATE-family multidrug resistance-associated protein, may aid in detoxification and antioxidant reactions. This process helps plants remove harmful substances and reduce oxidative damage. Their expression patterns at various low-temperature time points validate their significance in adaptation. Notably, their elevated expression under prolonged cold stress indicates pivotal roles in maintaining cellular stability and promoting repair. The strong correlation (0.78) between the Plum module and proline, along with enrichment in multiple metabolic pathways, reveals its multifaceted roles in low-temperature adaptation. Proline, a crucial osmotic regulator, accumulates under low-temperature stress, maintaining cellular osmotic balance and preventing dehydration and membrane damage. Additionally, proline acts as an antioxidant, directly scavenging ROS and mitigating oxidative damage [56]. KEGG enrichment of the Plum module showed significant enrichment in the alpha-linolenic acid metabolism, linoleic acid metabolism, MAPK signaling, glutathione metabolism, and peroxisome-related pathways. These enrichments suggest that the Plum module enhances plant adaptability to cold stress through synergistic metabolic and signaling interactions. Specifically, enrichment in alpha-linolenic and linoleic acid metabolism suggests the Plum module enhances membrane stability and functionality by adjusting fatty acid composition and fluidity, preventing temperature-induced membrane rupture [57]. MAPK signaling pathway enrichment underscores the module’s role in signal transduction under cold stress, suggesting plants regulate downstream genes through MAPK signaling to initiate stress resistance mechanisms, including antioxidant enzyme expression and cellular protection activation [58]. Glutathione metabolism and peroxisome-related pathway enrichment further highlight the Plum module’s role in antioxidant defense. By removing ROS, it reduces oxidative damage and enhances plant stress resistance [59]. The core gene EVM-20prediction-20ctg4298.8 and its interaction partners EVM-20prediction-20ctg1327.25 and EVM-20prediction-20ctg87.24 in the Plum module play central regulatory roles in cold stress responses. The core gene may enhance cold resistance by modulating MAPK signaling and fatty acid metabolism, maintaining membrane stability and effective signal transduction. Conversely, EVM-20prediction-20ctg1327.25 (MAPK signaling) and EVM-20prediction-20ctg87.24 (fatty acid metabolism) further support the Plum module’s multi-layered regulatory functions in cold adaptation. Heatmap analysis shows these core genes are rapidly upregulated after 2 h of cold exposure, indicating an early adaptive response. This indicates their key roles in initiating early cold adaptation. This dynamic pattern reveals the time-dependent regulation of Plum module genes in cold stress, highlighting their rapid response functions. The Pink and Plum modules work synergistically to adapt Caucasian clover to low temperatures, collaboratively regulating multiple critical biological processes.

The Pink module primarily protects cells and alleviates oxidative stress via secondary metabolite accumulation and antioxidant defenses, while the Plum module maintains membrane stability and intracellular osmotic balance through osmotic regulation and lipid metabolism. Together, these two modules enable Caucasian clover to effectively respond to physiological and metabolic stresses at low temperatures, ensuring cellular structural and functional integrity. Furthermore, core genes in both modules not only play pivotal roles in their own networks but also enhance cold resistance through cross-regulation and signal transduction. For example, EVM-20prediction-20ctg2798.293 in the Pink module may influence Plum module MAPK signaling by adjusting protein synthesis and signal transduction, thus facilitating a more efficient stress response. These inter-module interactions and synergistic effects illustrate how plants coordinate multiple biological processes via multi-level regulatory networks, achieving comprehensive adaptation to complex environmental stresses.

5. Conclusions

This study combines metabolomics, transcriptomics, and Weighted Gene Co-expression Network Analysis (WGCNA) to investigate the molecular mechanisms underlying Caucasian clover’s response to low-temperature stress. The results show that low-temperature stress regulates the cold tolerance of Caucasian clover via multiple biological pathways, primarily enhancing antioxidant defense, fatty acid metabolism, membrane lipid repair, and signal transduction. Metabolomic analysis revealed upregulation of fatty acid metabolism and amino acid synthesis, particularly the accumulation of antioxidant compounds, which are crucial for maintaining membrane stability and enhancing antioxidant capacity, thus protecting the plant from oxidative damage caused by low temperatures. Transcriptomic analysis confirmed these findings, showing significant upregulation of genes involved in antioxidant defense and membrane lipid repair, particularly those in the ABA signaling pathway and antioxidant enzymes, thus enhancing cold tolerance and supporting plant growth under low-temperature conditions. WGCNA identified several gene modules associated with low-temperature adaptation, particularly those involved in antioxidant defense, fatty acid metabolism, signal transduction, and membrane repair. The coordinated expression of these modules underscores their central role in cold stress adaptation. This study also investigated the isoflavonoid biosynthesis pathway, showing that low temperatures induced upregulation of key enzyme genes, such as Isoflavone Synthase (IFS) and Isoflavone-7-O-Glucosyltransferase (IF7GT), promoting antioxidant metabolite accumulation and further enhancing cold tolerance. This research significantly contributes to the understanding of the molecular mechanisms underlying plant adaptation to low-temperature stress, providing novel insights into the complex interactions between metabolic changes and gene expression that drive cold tolerance in plants. However, while these molecular insights are promising, evaluating the most effective approach for improving cold tolerance in breeding requires a careful, targeted analysis of the identified key genes and pathways. By selecting the most relevant genes and metabolic pathways, breeders can focus on genetic markers that have a direct impact on cold resistance, ensuring that future crop varieties are better equipped to thrive in cold climates.

In conclusion, this study offers novel molecular insights into the cold tolerance mechanisms of Caucasian clover, revealing the complex interplay between metabolite accumulation and gene expression. These findings contribute to the development of molecular markers and candidate genes for breeding cold-resistant crops. Future research should focus on validating these key genes in other cold-tolerant crops, advancing molecular breeding strategies for enhanced cold resistance and agricultural productivity.