Comparative transcriptomic and metabolomic analyses provide insights into the responses to high temperature stress in Alfalfa (Medicago sativa L.)

Abstract

High temperature stress is one of the most severe forms of abiotic stress in alfalfa. With the intensification of climate change, the frequency of high temperature stress will further increase in the future, which will bring challenges to the growth and development of alfalfa. Therefore, untargeted metabolomic and RNA-Seq profiling were implemented to unravel the possible alteration in alfalfa seedlings subjected to different temperature stress (25 ℃, 30 ℃, 35 ℃, 40 ℃) in this study. Results revealed that High temperature stress significantly altered some pivotal transcripts and metabolites. The number of differentially expressed genes (DEGs) markedly up and down-regulated was 1876 and 1524 in T30_vs_CK, 2, 815 and 2667 in T35_vs_CK, and 2115 and 2, 226 in T40_vs_CK, respectively. The number for significantly up-regulated and down-regulated differential metabolites was 173 and 73 in T30_vs_CK, 188 and 57 in T35_vs_CK, and 220 and 66 in T40_vs_CK, respectively. It is worth noting that metabolomics and transcriptomics co-analysis characterized enriched in plant hormone signal transduction (ko04705), glyoxylate and dicarboxylate metabolism (ko00630), from which some differentially expressed genes and differential metabolites participated. In particular, the content of hormone changed significantly under T40 stress, suggesting that maintaining normal hormone synthesis and metabolism may be an important way to improve the HTS tolerance of alfalfa. The qRT-PCR further showed that the expression pattern was similar to the expression abundance in the transcriptome. This study provides a practical and in-depth perspective from transcriptomics and metabolomics in investigating the effects conferred by temperature on plant growth and development, which provided the theoretical basis for breeding heat-resistant alfalfa.

Background

Temperature is one of the key environmental factors affecting the growth and development, and geographical distribution of plants, as well as quality and productivity. The entire plant life cycle is affected by environmental temperatures [1, 2]. Recently, global warming has led to frequent extreme high temperature stress events, and a global average surface temperature increase of 1.5 degrees Celsius to 4.8 degrees Celsius by 2100 [3], which have become one of the most important limiting factors to crop productivity and ultimately food security [4, 5], has led to crop yield losses, nutrition quality decreases, and widespread threats to food security and agricultural sustainability [6,7,8], raising a range of potential social problems. Elucidating the basic mechanism of plant response to high temperature, mining and identifying the gene loci of resistance to high temperature, and elucidating their functions and regulatory mechanisms are of great significance for solving the problem of food security caused by global warming.

High temperature stress not only affects the phenotype of plants, but also destroys cell homeostasis, seriously affecting the growth and development of plants, such as causes the reduction in net photosynthesis, leaf area, reduced biomass accumulation, and even contributing to their death [9, 10]. Research shows that three types responses have been documented in plants under HTS. HTS may alter metabolite event and inhibit various physiological and biochemical reactions of plants, including changing in water status, cell membrane stability, photosynthesis, secondary metabolites, plant hormone levels [11]. HTS can rapidly inhibit photosynthesis by changing the internal structure of the chloroplasts, reducing the abundance of photosynthetic pigments, and damaging photosystem II [12]. In addition, HTS can affect the protein synthesis and the stability of proteins and membrane, induce the accumulation of reactive oxygen species (ROS) [13], cause organelles to malfunction, alter plant hormone signaling [14], calcium signaling, lipid signaling and Kinases (MAKP, CBK, CDPKs), and induce transcriptomic reprogramming and metabolomic changes [15]. Therefore, it is important to unravel how plants respond to high temperature stress to improve the heat tolerance of plants.

In recent years, great progress has been made in the study of temperature stress, and the molecular regulatory network of plants to sense and respond to temperature stress has become increasingly rich. The molecular mechanisms of heat stress responses in plants are complex and controlled by multiple genes, genetic-based mechanisms of heat stress tolerance including the increased of stress signal, the control of transcription, the expression of heat shocked genes (HSG) and genes linked oxidative and osmotic stress [5].

A series of high temperature-related genes, small RNAs, and proteins have been reported. Plants have developed a set of sophisticated mechanisms that allow them to withstand HTS, the heat signal is perceived by putative heat sensors, such as the photoreceptor PHYB (photochrome B), ELF3 (EARLY FLOWERING 3) and PIF7 (photochrome interaction factor 7) [16]. PHYB is a temperature receptor of plant, which has been proved to play a certain role in regulating multiple processes such as light, temperature, hormones, stomatal opening and closing, and plant development [17], PHYB senses the heat signal by changing its protein activity (the active Pfr type changes to the inactive Pr type) [18, 19]. The EC (evening complex) composed of ELF3, ELF4 (EAERLY FLOWERING 4) and LUX (LUX ARRHYTHMO) is the core component of plant biological clock. EC not only regulates circadian clock gene expression, but also acts as an important environmental receptor to coordinate plant responses to the environment, was recently identified as a thermosensor [20]. The newest studies proved that ELF3 protein inhibits flowering gene expression, and ELF3 protein contains a protein-like domain (PrD) and PrD mediated ELF3 phase separation can sense ambient temperature changes, can unblock this transcriptional inhibition and promote flowering and growth, which is a new temperature sensing mechanism [21]. Importantly, Hsps (small heat shock proteins) bind unfolding proteins, playing a pivotal role in the maintenance of proteostasis in virtually all living organisms, which can help to prevent protein misfolding and aggregation at high temperatures and stabilize intracellular proteins and cellular homeostasis [22]. HSF (heat shock transcription factor), as one of the key components of plant heat stress response, plays an important role in plant heat stress memory. These HSPs are regulated by activating HSFs to enhance heat tolerance in plants. Members of the HSF family that function in heat stress memory include HsFA1, HsFA2, and HSFA3. Research finding HsFA1 function as ‘‘master regulators’’ by activating HSR genes and enhancing thermotolerance in Arabidopsis [23]. The negative regulatory component BIN2 of BR signaling pathway can phosphorylate and inhibit the nuclear localization and DNA binding ability of heat-shock transcriptional regulator HsfA1d, BR activates HsfA1d by disinhibiting the inhibitory effect of BIN2 on HsfA1d to enhance plant heat tolerance [24]. HsFA2 can be induced by HSFA1 protein and play a role specifically in heat stress memory, HSFA3 can effectively promote transcriptional memory by recruiting histone H3K4 methylation, either alone or in conjunction with HSFA2 [25]. In addition, studies found that various transcription factors and genes also play a role in the heat response mechanism, such as WRKY33/38/37/39, NAC019, SAG113, IAA29, PIF4/5, CBF2, BRI1, GT-1, HSP70, HSP90 [26,27,28]. These studies suggest that the regulatory mechanisms of plant tolerance to heat stress may involve multi-signaling and multi-regulatory gene.

Exploring the molecular mechanisms of plant response to heat stress can enhance the current understanding of genetic plant thermotolerance. While several regulatory pathways have been documented in heat stress response, there is little doubt that more regulatory pathways remain unexplored. A deeper comprehension of the stress response system could be achieved by integrating various approaches such as physiological, metabolomic, and transcriptomic data. The integration of transcriptomic and metabolomic data links the phenotype and physiology of plants, providing crucial information for further study on the intricate mechanisms of plant adaptation to high temperature stress. Zhao et al. [29] revealed that the thermotolerance of Paspalum wettsteinii might be related to the expression of heat-responsive genes enriched in energy metabolism, carbohydrate metabolism, heat shock protein, transcription factors, biosynthesis of secondary metabolites and antioxidant system, and the free fatty acids, amino acids, organic acids, flavonoids, and sugar accumulate under heat stress, respectively. Xie et al. [30] revealed that the regulation of the transcription factor family HSF and the purinergic pathway in response to high temperature stress to improve quinoa varieties with high temperature tolerance. The combined analysis of the transcriptome and metabolome can provide comprehensive information on the networks involved in regulating heat stress responses and provide a powerful analytical method for elucidating the changes from the genotype to the phenotype in plants under heat stress.

Alfalfa (Medicago sativa L.), a perennial legume forage that is widely cultivated for hay, pasture and silage productions. As a one of the most widely planted and economically valuable crops in the world, it has a high-quality forage due to its characteristics of high yield and excellent nutritional quality, which has the longest planting history and the largest planting area in the world [31, 32]. Temperature is one of the major environmental stressors that negatively affects alfalfa production. Global warming contributes to higher temperatures, significantly restricting alfalfa growth and production. At present, studies of high temperature stress on alfalfa have mainly focused on growth, physiological and biochemical responses. Although high temperature stress-related genes in alfalfa have been identified, there is no systematic study on the mechanisms underlying alfalfa response to heat stress. Based on the fact, to increase our understanding of response mechanisms to heat stress in alfalfa, a study on potted alfalfa plants was conducted in growing chambers under different temperature stress. Changes at the molecular level were investigated and linked to the adaptive responses of alfalfa including physiological, metabolic, and transcriptional traits. Understanding the physiological characteristics and identifying key genes and metabolites in alfalfa associated with high-temperature stress, providing the important reference for the study of high temperature stress in plants, and is vital for breeding new heat-resistant cultivars, especially with global climate change.

Materials and methods

Plant material and treatment

The alfalfa cultivars Zhongmu No.3 (ZM3) was used in the experiment, a high temperature-tolerant upland alfalfa genotype, it was assessed by Ningxia University for stress resistance. The healthy and plump alfalfa seeds were selected for disinfected the surface with 75% ethyl alcohol and germinate in an artificial climate chamber (BIC 250, Shanghai Boxun medical biological instruments CO., LTD, China) with a 14/10 h light /dark cycle at 25 ℃, then 6–8 germinated seeds were transplanted into plastic pots (140 mm×120 mm) containing a mixture of seeding substrate and vermiculite(4:1v/v)and grown in a greenhouse under14h /10 h light/dark at 70–75% relative humidity and average temperature of 25 ℃ [33]. Four-weeks-old seedlings were transferred to artificial climate chamber and treated with different temperature stress (25 ℃, 30 ℃, 35 ℃, 40 ℃) for 72 h. Each sample contained 8 seedings. After the treatments, the leaf of seedlings were harvested for transcriptomic (three independent treatment), metabolomic (four independent treatment) and physiological analyses (three independent treatment).

Physiological measurement of alfalfa

Catalase(CAT), peroxidase (POD), Ascorbate peroxidase (APx), Glutathione reductase (GR) and superoxide dismutase (SOD) activity, and the content of malondialdehyde (MDA), were determined according to the instructions of the reagent kit (Beijing Solarbio life sciences & technology Co., Itd, China), and the contents of six endogenous hormone auxin (IAA), gibberellin (GA), cytokinin (CTK), brassinosterol (BR), ethylene (Eth), and abscisic acid (ABA) were determined by the enzyme-linked immunosorbent assays (ELISAs) according the manufacture’s protocol (Shanghai Enzyme Linked Biotechnology Co., LTD, China). Three independent experiments were conducted [34].

RNA extraction and quantitative real-time PCR (qRT-PCR)

Total RNA was extracted by used EASYspin Plant RNA Kit (Aidlab Biotech, Beijing, China) and subsequently reverse transcription used Hifair^® III 1st Strand cDNA Synthesis SuperMix for qPCR (gDNA digester plus) from shanghai Yesean Biotech. 1000 ng RNA to reverse transcribe into cDNA, and diluted the cDNA four times for qRT-PCR. The qRT-PCR was performed by using the Gene Applied Biosystems^® 7500 Fast and TransStart Top Green qPCR SuperMix (TransGen Biotech, Beijing, China). The qRT-PCR was conducted in a 20 µL volume containing 1 µL of diluted cDNAs, 0.4 µL of each primer (showed in Table S9), 0.4 µL of passive reference dye, 7.8 µL H₂O, and 10 µL of Top Green qPCR SuperMix. Thermal cycling involved heating at 95 ℃ for 300 s, followed by 40 cycles of 95 ℃ for 15 s, 60 ℃ for 30 s and 72 ℃ for 30 s. The MsActin gene was used as an internal control, and the relative expression levels were calculated by the comparative 2^−△△Ct method [35].

Transcriptomic sequencing and data analysis

The cDNA library for transcriptomic sequencing was constructed and sequenced by APTBIO (Shanghai, China). Paired-end libraries were prepared using a ABclonal mRNA-seq Lib Prep Kit (ABclonal, China) following the manufacturer’s instructions. Adaptor-ligated cDNA was used for PCR amplification. PCR products were purified (AMPure XP system) and library quality was assessed on an Agilent Bioanalyzer 4150 system, Then the cDNA library was sequenced on the DNBSEQ-RST7-T7 high-throughput sequencing platform, and 150 bp raw FASTQ data were obtained. Quality filtering was performed using Fastp [36], and high-quality reads were then mapped to the genome using HISAT2 [37]. Differential expression analysis was performed using the R package DESeq2 [38]. Genes with the standard of a |Log2FoldChange| ≥ 1 and Padj < 0.05 were identified as differentially expressed genes (DEGs) [39]. The BLAST2GO and Kobas were used for the annotations of GO and KEGGs, respectively, and cluster Profiler R software package was used for enrichment analysis [40]. The R and TBtools were used for data visualization [41].

Sample extraction and measurements for metabolomic analysis

The frozen samples were ground in liquid nitrogen, and aliquots of 80 mg weight for metabolite extraction, then adding 1000 µL methanol/acetonitrile /water (2:2:1, v/v/v), followed by vorteing for 30 s, low temperature ultrasonic 30 min, and preserved at -20 ℃ for 10 min [42]. The mixture was centrifuged for 20 min (14,000 g, 4 °C), and then dried in a vacuum centrifuge. The samples were redissolved in 100 µL of acetonitrile/water (1:1, v/v) solvent and centrifuged for 15 min (14,000 g, 4 ℃), and the supernatant was collected and filtered for UPLC–MS/MS analysis. The metabolites were detected by APTBIO (Shanghai, China) and analyzed using a UHPLC (Vanquish UHPLC, Thermo) coupled to an Orbitrap.

Differential metabolite analysis

After sum-normalization, the R package were used to analyze the processed data. Multivariate analysis of the data was performed, including Pareto-scaled principal component analysis (PCA) and orthogonal partial least-squares discriminant analysis (OPLS-DA). Student’s T-test and variable influence on projection (VIP) values were utilized to determine the statistical significance of identified metabolites. Significant changes in metabolites were screened based on OPLS-DA results, with VIP > 1 and a p value < 0.05. Pearson’s correlation analysis was conducted to assess the relationship between variables. Significantly regulated metabolites between groups were determined by VIP > 1, p value < 0.05, and fold change (FC) > 1.2 or < 0.83 [43]. Identified metabolites were annotated and mapped to KEGG pathway database (http://www.kegg.jp/kegg/pathway.html). Pathways with significantly regulated metabolites mapped to MSEA (metabolite sets enrichment analysis), and their significance was determined by the p-values of hypergeometric test [44].

Integrated analysis of transcriptome and metabolome

The orthogonal partial least squares discriminant analysis was complete by the APTBIO. Query all DEGs and DEMs from the online KEGG (http://www.kegg.jp/) and mapped to pathways to get the common pathway information. The KEGG annotation and enrichment results for the two omics were combined using the statistical program R (version 3.5.1). The screening criterion for significant enriched pathway were set with a P-value < 0.05.

Statistical analysis

The Origin 2023 was used to plot charts. Significant differences among all the experimental treatments were determined using a One-way analysis of variance (ANOVA) followed by Tukey’s test by IBM SPSS version 26.0 software. The free online platform APTBIO Cloud Platform, Omicshare tools (https://www.omicshare.com/tools/), and Major Platform ( http://cloud.majorbio.com/page/tools.html ) were used to draw heat maps, GO enrichment scatter and KEGG enrichment scatter. The merge images were used Adobe Illustrator 2022.

Results

Physiological characteristics of under temperature stress

HTS disrupts the balance of ROS scavenging and lead to oxidative stress. Therefore, physiological characteristics of Medicago sativa L. under three kinds of temperature stresses were determined. The results showed that antioxidant enzyme activities changed and the balance of ROS scavenging was destroyed under all three stresses. The activities of GR (Fig. 1b), SOD (Fig. 1c), POD (Fig. 1d), CAT (Fig. 1e), APX (Fig. 1f) were significantly increased under T30, T35, T40, and the CAT activity increased higher than other enzymes (Fig. 1e). The results of The accumulation MDA decreased under all three stresses (Fig. 1a), was negatively correlated with enzyme activity. Endogenous hormone plays a significant role in plant stress response, this study determined six endogenous hormones. The results showed that the endogenous hormone fluctuated under different treatments, the content of CTK (Fig. 1g), IAA (Fig. 1h), ABA (Fig. 1i), GA (Fig. 1k) was higher than the control significantly, especially in T35 and T40 treatment. ETH (Fig. 1j) increased slowly under three treatments, but till lower than CK. With the temperature gradient, the content of BRs (Fig. 1l) increased, decreased and then increased, and the changes were significantly higher than CK under T40 treatments. The above evidence indicated that the responses to three temperature in alfalfa, the membrane damage caused by T40 stress is more serious, while the activity of enzyme and content of endogenous hormone were highest than other, this can effectively clear up reactive oxygen species, and alfalfa was protected to better adapt to temperature stress.

Fig. 1

Physiological characteristics of alfalfa under different temperature stress. (a) MDA content; (b–f) GR (b), SOD (c), POD (d), CAT (e), APx (f) activities; (g–l) CTK (g), IAA (h), ABA (i), ETH (j), GA (k), BR (l). a, b, c, d represents a significant difference (P < 0.05)

Transcriptomic analysis of alfalfa during high temperature stress

In order to understand the gene expression of alfalfa in responses to temperature stress, we performed RNA-seq analysis. High quality reads were obtained. After quality control filtration, the average Q20 was 98.23%, and the average Q30 was 95.11%. The mean GC content was 44.88% (Table S1), indicating that the obtained data were useful for further analysis. The trends of separation obtained from the principal component analysis (PCA) of the samples reflected the degree of differences among them. The PCA results of the RNA-seq samples showed that PC2 separated the control group (CK) from the treatment groups (T30, T35 and T40), and PC1 separated the treatment groups (Fig. 2a). DEGs were identified under T30, T35 and T40 stress to explore the similarities and differences at the transcription level in response to the three stresses in alfalfa. In total, 3, 418, 5, 483 and 4, 341 DEGs were identified in alfalfa during T30, T35 and T40 stress, respectively (Fig. S1). Among them, 1, 876, 2, 815 and 2, 115 genes were up-regulated, while 1, 542, 2, 667 and 2, 226 genes were down-regulated (Fig. 2c). In addition, 1, 810, 2, 889 and 2, 410 genes were specifically differentially expressed under T30, T35 and T40 treatment, respectively, while 544 genes were differentially expressed after all treatments (Fig. 2b). Collectively, these data showed that more genes were up-regulated than down-regulated in response to T30 and T40 stress, more down-regulated genes in response to T35 stress.

Fig. 2

RNA-seq analysis of alfalfa during different temperature stress. (a) PCA of RNA-seq samples. (b) UpSet plot of DEGs of alfalfa under temperature stress. (c) Statistical analysis of DEGs; up represents upregulated genes, and down represents downregulated genes

To further explore the DEGs involved in the responses to different temperature stress, Gene Ontology (GO) enrichment analysis was performed. All (up and down genes) DEGs in T30_ vs_ CK, T35_ vs. _CK, T40_ vs_ CK treatment groups were categorized into three GO classification (including 228 subclassification) according their functional annotations, specifically, biological process (BP), molecular function (MF) and cellular component (CC) (Fig. 3). In addition, most DEGs appeared in the classification of BP and MF groups, and most DEGs belong to MF group gathered into “binding” and “catalytic activity” subgroups. Under T30, T35 and T40 stress, DEGs significantly enriched in the biosynthetic process, organic substance biosynthetic process, and metabolic process.

Fig. 3

GO enrichment analysis of DEGs that were specifically expressed in response to different temperature stress. (a–c) T30, T35 and T40 stress. The dot area represents the number of genes, and the color represents Padjust value

KEGG enrichment analysis of EDGs in three comparison groups was performed, there were 59, 102, and 80 KEGG pathways identified in T30 vs. CK, T35_ vs. _CK, and T40_ vs. _CK comparison group. The pathway of oxidative phosphorylation, glycerolipid metabolism, thermogenesis, glycerophospholipid metabolism, tropane, piperidine and pyridine alkaloid biosynthesis were significantly enriched under T30_ vs. _CK (Fig. 4a), the pathway of carbon fixation in photosynthetic organisms, glyoxylate and dicarboxylate metabolism was significantly enriched under T35_ vs. _CK (Fig. 4b), the pathway of plant hormone signal transduction, carbon fixation in photosynthetic organisms, glyoxylate and dicarboxylate metabolism, photosynthesis-antenna proteins, ascorbate and aldarate metabolism were significantly enriched under T35_ vs_ CK (Fig. 4c).

Fig. 4

KEGG enrichment analysis of DEGs in response to different temperature. (a–c) T30, T35 and T40 stress. The dot area represents the number of genes, and the color represents Padjust value

To analyze the regulatory mechanisms of high temperature responsive genes, 86 (47 up-regulated and 39 down-regulated), 142 (92 up-regulated and 50 downregulated) and 131 (65 up-regulated and 66 down-regulated) DEGs encoding TFs were identified in Medicago Sativa L. (Table S2, S3, S4). A total of 19, 24, and 19 transcription factor families were identified in T30,T35, and T40, including AP2, zf-Dof, HSF_DNA bind, CSD, Zim, Homeobox and WRKY transcription factor families (Fig. 5a, b, c), The major transcription factor families three were AP2 (range from 19.9 to 46.4%), zf-Dof (range from 14.4 to 20.9%), HSF_DNA bind (range from 2.3 to 8.1%). The results suggest that these transcription factors may be the main high temperature-responsive genes and play important roles in alfalfa high temperature tolerance.

Fig. 5

Distribution of different temperature stress responsive transcription factor genes. (a–c) Classification of TFs under T30, T35, T40 stress, and the TFs that accounted for less than 1% were summarized to others

Metabolomic analysis of alfalfa during high temperature stress

The positive (POS) and negative (NEG) ion modes were used for qualitative and quantitative evaluation to improve the coverage of metabolites and detection effect. The PCA results of metabolic samples were consistent with the RNAseq results, which revealed high similarity among the four biological replicates under each treatment. There were differences in the samples analyzed in the three cases: PC1 separated the T40 treatment t from the T30, T35 and CK samples, and PC2 separated the CK sample from the T30 and T35 samples (Fig. 6a). Under comprehensively unbiased non-targeted metabolomic profiling, a total of 2, 121 metabolites were identified by positive (1, 223) and negative (898) ion patterns, and grouped into 16 classes (Table S5). According to the metabolite classification results, there were 7 types of compounds that account for more than 80% of the total metabolites in leaves, including lipids and lipid-like molecules (22.68%), organic acids and derivatives (20.23%), organoheterocyclic compounds (13.20%), benzenoids (12.82%), Phenylpropanoids and polyketides (7.54%), organic oxygen compounds (7.26%), organic nitrogen compounds (1.84%), nucleosides, nucleotides and analogues (1.51%).

The overall distribution of differential metabolites in each comparison group is shown in Fig. 6c. There were 247 differential metabolites in T30_vs_CK comparing group, of which 174 were significantly up-regulated and 73 were significantly down-regulated, respectively. T35_vs_CK comparison group had 245 differential metabolites, among which 188 and 57 were significantly up-regulated and down-regulated, respectively. There were 250 differential metabolites in the T40_vs_CK comparing group (197 significantly up-regulated and 53 down-regulated).

Among these metabolites, 76 (49 up-regulated and 27 down-regulated) metabolites were differentially expressed under three temperature stress, and 97, 71 and 126 DEMs were specifically expressed under T30, T35 and T40 stress, respectively (Fig. 6b). In general, the number of up-regulated metabolites was higher than that of down-regulated metabolites, especially under T40 stress. The DEMs identified during including lipids and lipid-like molecules, benzenoids, organic acids and derivatives, organoheterocyclic compounds, organic oxygen compounds, organic nitrogen compounds, Nucleosides, nucleotides and analogues.

The KEGG enrichment analysis revealed that differential metabolites identified in T30 _ vs. _CK, T35_vs_CK, and T40_vs_CK treatment groups significantly enriched in 11, 23, and 16 pathways related to plant metabolism (Tables S6, S7, S8). The significantly enriched KEGG pathways for differential metabolites shared by three treatment groups include ABC transporters; plant hormone signal transduction; alanine, aspartate and glutamate metabolism; glyoxylate and dicarboxylate metabolism; citrate cycle (TCA cycle); glycine, serine and threonine metabolism. Consequently, temperature stress significantly changed some metabolites generation in alfalfa leaves and markedly modified the enrichment route of these metabolites. In addition, we found that another 9 pathways, including protein digestion and absorption; mineral absorption, pantothenate and CoA biosynthesis, carbon metabolism, cysteine and methionine metabolism, galactose metabolism, cAMP signaling pathway, beta-Alanine metabolism, and phenylpropanoid biosynthesis were all noteworthily altered after the leaf was exposed to different temperature. Therefore, temperature treatments influence the biosynthesis of some secondary metabolites, such as sucrose, betaine, biotin, 2-cis-4-trans-abscisic acid, salicylic acid, trans-zeatin, glyceric acid, isocitrate, succinate, glutamic acid etc. In conclusion, these bioactive substances and changes in their enrichment pattern suggested that alfalfa adapts to temperature stress by activating the biosynthesis of some secondary metabolites.

Fig. 6

Metabolic analysis of alfalfa under different temperature stress. (a) PCA of metabolic samples. (b) Venn diagram showing the shared and unique DEGs between T30, T35 and T40 stress. (c) Statistical analysis of DEMs. the up-regulated metabolites and down-regulated metabolites were shown in the color of orange and green, respectively, and the numbers in the bar represent the number of metabolites

Integrated analysis of transcriptomic and metabolomic expression levels

Integrated metabolomic and transcriptomic methodology was used to assess the potential associations between DEGs and DEMs to reveal the pathways through which they are jointly involved. There were 30, 33 and 34 transcripts and metabolites common shared under T30, T35, T40 treatments (Fig. 7a-c). The top 10 KEGG pathways with the most participation of gene and metabolite identified in this experiment are illustrated in Fig. 4, and the results showed that the KEGG pathways in which differential metabolites and genes were jointly involved were glycerolipid metabolism (ko00561), glycerophospholipid metabolism, Glyoxylate and dicarboxylate metabolism (ko00630), Plant hormone signal transduction (ko04705), respectively (Fig. 7). Several DEGs and the differential metabolites are involved in all three pathways, although the DEGs did not directly regulate the production of these different metabolites.

Furthermore, the results indicated that the DEGs were significantly enriched in oxidative phosphorylation and glycerolipid metabolism, while the DEMs were significantly enriched in the glyoxylate and dicarboxylate metabolism and plant hormone signal transduction pathway during T30 stress (Fig. 7d). Under T35 stress, DEGs were significantly enriched in the glyoxylate and dicarboxylate metabolism and glycerolipid metabolism pathway, and the DEMs were enriched in the citrate cycle (TCA cycle) and plant hormone signal transduction pathway (Fig. 4e). During T40 stress, both DEGs and DEMs were significantly enriched in the plant hormone signal transduction and glyoxylate and dicarboxylate metabolism pathway (Fig. 7f).

Fig. 7

Correlation analysis between metabolomics and RNA-seq data. (a–c) venn diagram of pathways involving differential genes and differential metabolites under T30, T35, T40 stress; (d–f) The top 10 pathways contain the most genes and metabolites involved together under T30, T35, T40 stress, the x-axis indicates the shared KEGG pathways, and the y-axis indicates the number of genes and metabolites involved in the KEGG pathways

Analysis of shared or specific DEGs and DEMs under different temperature stress

According to the integrated analysis of transcriptomic and metabolomic, the top 10 KEGG were analysed under T30, T35 and T40 stress. A total of 74 DEGs were shared in top 10 KEGG (Fig. 8), among which 32 and 42 DEGs were up-regulated and down-regulated under T30, 34 DEGs up-regulated and 40 down-regulated under T35, 15 DEGs up-regulated 59 down-regulated under T40. The 74 DEGs involved 16 KEGG pathway, among which 28 DEGs belongs to plant hormone signal transduction (ko04075), 15 DEGs belongs to glyoxylate and dicarboxylate metabolism (ko00630), 10 DEGs belongs to glycerolipid metabolism (ko00561), 9 DEGs belongs to carbon fixation in photosynthetic organisms (ko00710). A total of 41 metabolites were differently expression under temperature stress in top 10 KEGG (Fig. 9a), which 10 metabolites shared under the T30, T35 and T40 stress treatments (Fig. 9b).

In addition, KEGG results of the DEMs and DEGs showed that DEMs were associated with glyoxylate and dicarboxylate metabolism (ko00630), plant hormone signal transduction (ko04075), arginine and proline metabolism (ko00330), glycine, serine, and threonine metabolites (ko00260), cysteine and methionine metabolism (ko00270). Interestingly, among the pathways related to DEGs, the plant hormone signal transduction and glyoxylate and dicarboxylate metabolism were also shared with the DEMs, but the other pathways were different. Specifically, DEGs and DEMs were also related to phenylalanine metabolism, oxidative phosphorylation, and ascorbate and aldarate metabolism. The results demonstrated that under high temperature stress, DEGs participated in the synthesis of stress-related secondary metabolites, which may reduce stress damage through the antioxidant pathway.

The 18 specific DEMs indicated differences among the responses of alfalfa to different temperature stress (Fig. 9c). A total of 8 metabolites (abscisic acid, salicylic acid, L-hydroxyarginine, phosphocreatine, .alpha.-keto-.gamma.-(methylthio) butyric acid1, betaine, phenylalanine, dexpanthenol) specifically differentially expressed involed 6 KEGG pathway during under T30 (Fig. 9c), including the arginine and proline metabolism (Ko00330), pantothenate and CoA biosynthesis (ko00770), protein digestion and absorption(ko04974). There were 2 DEMs (glycerophosphocholine, Acetylcholine) specifically expressed involved in glycerophospholipid metabolism (ko00564) under T35 stress. There were 8 DEMs (Jasmonic acid, glutamic acid, D-proline, glycine, L-aspartic acid, L-threonine, L-methionine, D-galacturonic acid) were involved in plant signal transduction (ko04075), glyoxylate and dicarboxylate metabolism (ko00630), cysteine and methionine metabolism (ko00270), ascorbate and aldarate metabolism (ko00053), arginine and proline metabolism (ko00330) during T40 stress. These specific DEGs were involved in the response to high temperature stress by participating in various pathways to induce the synthesis of specific metabolites.

Fig. 8

Analysis of differentially expressed genes involved in KEGG pathway under T30, T35 and T40 stress(Log₂FC was used to construct the heatmap)

Fig. 9

Analysis of differentially expressed metabolites involved in the response to T30, T35 and T40 stress(Log₂FC was used to construct the heatmap). (a) Heatmap analysis showing the expression trend of DEMs; (b) Analysis of DEMs shared expressed duringT30, T35 and T40 stress; (c) Analysis of DEMs specifically expressed during T30, T35 and T40 stress

Plant hormone signal transduction

Zeathin biosythesis, caroterniod biosythesis, Phenylpropanoid biosythesis, α-linolenic acid biosynthesis and brassionoteriod biosynthesis were significant changed in the plant hormone signal transduction pathway (Fig. 10a). The metabolite cytoskinine were significantly upregulated in zeanthin biosythesis, and the gene CRE1 (MS. gene68824, MS. gene65393) was downregulated. The metabolite jasmoruic acid were significantly upregulated in the α-linolenic acid biosynthesis, with the gene JAR1 (MS. gene22809) and COI1 were upregulated, and JAZ (MS. gene89156; MS.gene09703; MS.gene54190; MS.gene21515; MS.gene26683) was downregulated (Fig. 10b). The metabolite aslicylic acid were significantly upregulated in the Phenylpropanoid biosynthesis, with gene TGA (MS. gene29028) was downregulated (Fig. 10b). In addition, the gene AUX1 (MS.gene 90982; MS.gene33487) was down-regulated in Tryptophan biosythesis, the gene PYR/PYL4 (MS.gene33704; MS.gene33451; MS.gene67442; MS.gene95667) and BZR1/2(MS.gene59331) were up-regulated (Fig. 10b).

Fig. 10

Simplified representation of DEGs and DEMs aboundance abundance in hormone signal transduction metabolism pathway. (a) Simplified representation of hormone signal transduction metabolism pathway. Transcripts (substances on either side of the arrow) and metabolites(circle) with different abundance involved in plant hormone signal transduction were mapped to corresponding metabolic pathways of KEGG. Red indicates up-regulation, and green indicates down-regulation, blue indicated unchanged. (b–c) The genes and metabolites implicated in plant signal transduction pathway, and Log₂FC values of each gene and metabolites

Glyoxylate and dicarboxylate metabolism

In the metabolism pathway of glyoxylate and dicarboxylate, the citrate, Cis-aconitate (EINECS), I-socitrate, succorate and L-Glutanate were significantly upregulated, while D-glycerate, glutamine synthetase (GS), ribulose-bisphotosphate carboxylase large chain (RuBisCO), D-glycerate 3-kinase (GLYK) were remarkably downregulated (Fig. 11a). The identified genes involved in encoding EINECS, GS, RubisCO, and GLKY are shown in the Fig. 11. The significant upregulation of EINECS was beneficial to the increment of I-socitrate and succorate, by contrast, the conspicuously down-regulated GS was not conducive to the generation of L-glutamine (Fig. 11a). Over all, the evident down-regulation of RuBisCO and GLYK was disadvantageous to getting D-glycerate and 3-phospho-D-glycerate, respectively.

Fig. 11

Simplified representation of DEGs and DEMs abundance in glyoxylate and dicarboxylate metabolism pathway. (a) Simplified representation of glyoxylate and dicarboxylate metabolism pathway. Transcripts (substances on either side of the arrow) and metabolites(circle) with different abundance involved in glyoxylate and dicarboxylate metabolism pathway were mapped to corresponding metabolic pathways of KEGG. Red and green means up-regulation and down-regulation, blue indicated unchanged. (b–c) The genes and metabolites implicated in glyoxylate and dicarboxylate metabolism pathway, and Log₂FC values of each gene and metabolites

qRT-PCR verification of RNA-seq data

To validate the RNA-seq results, randomly selected 9 DEGs for qPCR analysis. The expression levels of MS.gene08315, MS.gene59331, MS.gene89156, MS.gene73493 increased at T30 (relative expression level > 1) but decreased at T35and T40 (relative expression level < 1). Furthermore, although the expression levels of MS.gene54468, MS.gene92218, MS.gene54467, MS.gene80967, MS.gene76907 decreased (relative expression level < 1) with temperature change (T30 or T40), the expression levels increased at T35. The results showed that the expression trends of these genes obtained by the two technologies were basically consistent, indicating the validity of the RNA-seq results (Fig. 12).

Fig. 12

Comparison between the results of qRT-PCR and RNA-seq of selected genes. The blue bar represents FPKM value in RNA-seq data, the orange bar represents the qRT–PCR results and was calculated using the 2^−ΔΔCt method

Discussion

Temperature is a critical factor that significantly impacts agricultural development. In recent years, crop failures have been observed worldwide due to elevated high temperatures. HTS can lead to changes in plant morphology, physiology, and gene expression, ultimately inhibiting the growth and development of plant [45]. When plants are exposed to high temperature stress, they undergo a series of metabolic alterations and produce more metabolites to adapt to the ever-changing environment. These adaptations include the accumulation of osmotic adjustment substances (such as soluble protein, soluble sugars, proline, and betaine ) [46], changes in hormone levels (including ABA, IAA, CTK, GA, Eth, BRs, JA, SA ) [47] and alterations in antioxidant enzyme activity (such as SOD, POD, CAT, APx, GR) [48]. The metabolic regulation of alfalfa leaves was characterized by significantly up-regulated and down-regulated metabolites (such as jasmonic acid, abscisic acid, salicylic acid, citrate, isocitrate, succinate, glutamic acid) and alteration in some important metabolic pathways (a-linoleruic acid, tryptophan, phenylalanine, zeathin, caroterinoid, brassinoteriod) after exposure to 30 ℃, 35 ℃ and 40 ℃. These altered metabolites act as stress coordinators, alleviating the changes of alfalfa, and ultimately, reestablishing the homeostasis between the plant and its surroundings. At the same time, these genes/transcripts play an important role in regulation of plant physiology as precursors for the biosynthesis of various enzymes and metabolites.

Transcriptome analysis showed that a total of 13,242 DEGs were identified under HTS, with 6,806 up-regulated and 6,436 down-regulated DEGs. Among these identified DEGs, most DEGs were annotated to the plant hormone signal pathway, and most of DEGs were down-regulated under HTS; At the same time, the content of 5 metabolites mapped to plant hormone signal pathway accumulated or decreased under high temperature treatment. These results indicated that plant hormone signal pathway was regulated by high temperature in alfalfa. Based on the transcriptomic and metabiotic analysis, a lot of transcripts annotated to plant hormone signal transduction (ko04705), glyoxylate and dicarboxylate metabolism (ko00630), glycerolipid metabolism (ko00561), glycerophospholipid metabolism under high temperature stress. Consistent with the transcriptome data, qRT-PCR demonstrated five transcripts were less expressed, and four transcript was more expressed under high temperature stress in these three biosynthesis pathway among these transcripts. The qPCR results confirmed the changes in the expression levels of nine genes under high temperature stress.

JA is one of the important plant hormones, that plays an important role in regulating plant biological and abiotic stress response and growth and development [49]. Studies have shown that JA was a key hormone that tolerate bright light and high temperatures, respectively [50]. It is an oxidized fatty acid derivative in plants, synthesized mainly from α-linolenic acid via the octadecanoic acid pathway. The plant hormones JA and methyl jasmonate (MeJA) are vital signal transducers. Research conducted by zhang et al. [51] showed the tea HSFA1b/A2 can bind the JA inhibitor CsJAZ6 gene promoter and negatively regulate flavonoid synthesis in response to high temperature stress; Tian et al. [52] found the expression of JA synthesized gene TaOPR3 was induced by high temperature, over expression of TaOPR3 can significantly improve the high temperature tolerance of wheat, while reducing expression can reduce the high temperature tolerance of transgenic wheat. In summary, the JA metabolite in alfalfa leaves increased significantly after HTS in the α-linolenic acid metabolism pathway. In contrast, despite the down-regulation of JAR1 and COI of JAZ transcripts, stress-responsive signaling transmitted by jasmonate-related metabolites unchanged, and JAZ transcripts up-regulated. There seems to be a deep fine-tuning mechanism to support the trade-offs between growth and defense. Therefore, it is necessary to explore the deep regulation mechanism of this metabolic pathway.

SA is the central defense hormone in plants [53], through Isochorismate synthase (ICS) pathway and phenylalanine ammonia lyase (PAL) pathway biosynthesis [54]. SA’s main roles in plant tolerance to stresses can be summarized as following: the activation of plant defense mechanisms, the improvement of photosynthesis, the regulation of plant normal growth and development, and the enhancement of plant tolerance to drought, cold, heat, salt, and heavy metal stresses [55]. For all plants, temperature sensitivity appears to be specific to the salicylic acid signaling pathway, and high temperatures down-regulate the expression of salicylic acid response genes. Increasing temperature can inhibit the expression of ICS1 (a key gene for SA biosynthesis), and reduce the accumulation of salicylic acid at high temperature [56]. Previous research indicated NONEXPRESSOR OF PATHOGENESIS-RELATED GENES (NPRs) and TGACG mode-binding transcription factors (TGAs) mediate SA signaling [57], TGA1 and TGA4 activate the expression of SARD1 gene and indirectly regulate the biosynthesis of SA [58], and BR signals coordinate SA responses mediated by TGAs [59], respectively. In this study, the SA metabolite in alfalfa leaves increased significantly under HTS with TGA translator up-regulated, the gene’s function needs to be depth study.

As a key hormone in plant response to stress, ABA not only has the physiological function of activating NADPH oxidase (Rboh) and other enzymes to produce ROS, but also is closely related to sugar transport and starch anabolism in plant reservoir/source organs [60]. ABA and receptor PYR/PYL/RCAR protein constitute stable complex PYR/PYL/RCAR, binding and inhibited negative regulator 2 C protein phosphatase (PP2C), and released positive regulator SNF1-related protein kinase 2 (SnRK2), and the activated SnRK2 phosphorylates downstream transcription factors and ion channels, activating ABA signaling pathways and stress response processes [61, 62]. The gene PYR/PLY up-regulated in carotenoid biosynthesis in this study, it is speculated that PYR/PLY induced by HTS and positive regulated the ABA-mediated response to temperature stress in Medicago sativa L., but the regulatory mechanism of PYR/PYL can be further studied in the future.

BRs are a class of sterol hormones that were the first hormones identified as essential for plant growth and development, which has been widely studied in recent years [63]. Subsequent studies gradually found BRs has multiple functions, including the regulation of plant photomorphogenesis, root and stomatal development, abiotic stress and defense response [64]. It has been confirmed that BRs plays a positive regulatory role in response to heat stress [65]. Tae-Woo et al. [66] found BR signaling pathway is involved in ABA-mediated stomatal closure through CDL1, BRI1 and BAK1, and BAK1 can interact with OST1. As a kind of regulating various growth and development processes of plant hormone, is the regulation of stomatal movement by BR a strategy to balance plant growth and development and energy metabolism? More research is needed to answer these questions.

High temperature stress can cause the explosion of ROS in plants, resulting in irreversible cellular oxidative damage. Effective detoxification of ROS is believed to play a key role in enhancing plant tolerance to abiotic stresses [67]. Under high temperature stress, the plants will start the antioxidant system and produce a series of antioxidant enzymes, such as SOD, APX, POD, CAT, GR, AsA, and DHAR to clear the damage of ROS to the cell membrane and improve the heat resistance of plants. Our transcriptome analysis has identified genes and metabolic pathways related to the defence system responses to heat stress, and also found that the activity of antioxidant enzymes gradually increased with the increase of temperature stress, which also confirmed that the antioxidant oxidase system is one of the mechanisms of alfalfa in response to high temperature stress. SOD, CAT, POD, GR and APX protein acts as a crucial antioxidant enzyme in ROS metabolism in plants response to biotic and abiotic stresses, it can catalyze the transformation of H₂O₂ to O₂ and H₂O, and then remits the damage of H₂O₂ to plant cells. Accumulating evidences have indicated that APX and SOD to be involved in various abiotic stress related responses such as heat, drought, salt, low temperature conditions [68]. Yan et al. [69] found overexpression of CuZnSOD and APX enhance salt stress tolerance in sweet potato. Zhai et al. found overexpression of SIGR in transgenic tobacco can decrease the accumulation of ROS and increase the expression and activity of GR, APX, SOD and CAT under salt stress [70]. Shen et al. [71] found the instantaneous expression of SmAPX2 in vivo alleviates the accumulation of H₂O₂ under high temperature stress, the silencing of SmAPX2 decreased the tolerance of eggplant to high temperature stress, and was accompanied by downregulation of the transcription of genes related to high temperature defense and aba response, decreased APX enzyme activity, and increased H₂O₂ content. These studies indicated that antioxidant enzyme genes function in the process of growth, development, and response to abiotic stresses in plants. However, the functions and mechanisms of SOD, CAT, POD, GR, APX genes in alfalfa response to high temperature stress are largely indistinct. The mechanism of synergistic response of antioxidant enzymes and hormones to high temperature stress, and how plants respond to stress based on gene regulatory network need to be further studied in the future.

Conclusion

In the present work, we sequenced and assayed leaf transcriptome sequence and metabolome in alfalfa subjected to different temperature (30 ℃, 35 ℃ and 40 ℃) stress. The integrated metabolomics and transcriptomics analysis confirmed that the expression of some genes and the generation of some metabolites in alfalfa leaves have changed under stress relative to the control (CK). Additionally, the changes that occurred in two significantly enriched KEGG pathways (viz., plant signal transduction, glyoxylate and dicarboxylate metabolism) were analyzed through the integrated metabolomics and transcriptomics, where in some differential genes and metabolites participated in the pathways. Our results is expected provide valuable information for explore the high temperature regulation mechanism of alfalfa. However, the deep and comprehensive mechanism of genes and metabolites variation induced by high temperature stress is deserved further exploration. Further, coordinated multi-layer efforts including physiological, chemical, gene function, transgene, and approaches based on omics such as proteomics, scRNA-seq, lipidomics and nanobiotechnology are critical, and may yield more desirable results that reveal the effects of high temperature stress on the growth of plants from a new perspective in the future.