Abstract

   In this study, we aimed to investigate the heat tolerance mechanism in
   Rhododendron henanense subsp. lingbaoense (Rhl). Rhl seedlings were
   treated at 40℃ (RLH), 32℃ (RLM), and 24℃ (RLC), and the changes in
   transcriptome and metabolome were compared. Overall, 78 differentially
   expressed metabolites were detected, and 8450 differentially expressed
   genes (DEGs) were identified. KEGG analysis revealed that the DEGs in
   RLH vs. RLC were mainly enriched in photosynthesis, secondary metabolic
   biosynthesis, and flavonoid biosynthesis. Most genes encoding
   glutathione-S-transferase were upregulated, whereas genes related to
   heat shock proteins were significantly downregulated. 31 genes related
   to photosynthesis were significantly upregulated (P-value < 0.001). It
   was speculated that these DEGs are related to the response of Rhl to
   high temperature stress (HTS). Overall, 9 TF families might be the key
   regulators of Heat stress response pathways in Rhl. Mining of DEGs
   revealed that the expression of some genes related to heat stress
   function increased highly significantly, e.g., the Rhe008987 related to
   Glutathione-S-transferase, Rhe016769 encoding peroxidase, and Rhe001827
   encoding chalcone and stilbene synthases. Metabolome and transcriptome
   correlation analysis revealed that three comparison groups (RLH vs.
   RLC, RLH vs. RLM, and RLM vs. RLC) shared 12 metabolic pathways in
   which the DEMs were enriched. HTS inhibited or induced expression of
   genes in flavonoid biosynthesis pathway and led to decreace in
   kaempferol content and quercetin accumulation. HT induced expression of
   genes in ABC pathway, which may be one of the reasons for the
   significant accumulation of L-isoleucine, L-leucine, and L-proline. In
   this study, DEGs mining found that the expression of some genes related
   to heat stress function increased highly significantly. And two omics
   correlation analysis revealed that 12 metabolic pathways were enriched
   in three comparison groups. These results helped in elucidating the
   molecular mechanisms of response of Rhl to HTS.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s12870-025-06305-3.

   Keywords: Rhododendron henanense subsp. lingbaoense, High-temperature
   stress, Transcriptome, Metabolome, Flavonoid biosynthesis,
   Transcriptome sequencing

Introduction

   Rhododendron is a general name of plants from Rhododendron L. genus
   from Ericaceae family. These ornamental plants are highly valued
   worldwide [[42]1]. China is rich in wild rhododendron resources, with
   more than 600 species [[43]2], However, China’s rich rhododendron
   resources have not been optimally used in urban landscape construction
   and economic development. One of the important reasons for the hindered
   development of rhododendron industry in China is poor heat resistance
   of rhododendron [[44]3]. With the continuous social and economic
   developments, some wild plant species with small populations and little
   attention have been extinct or endangered due to the serious
   disturbance and narrow habitat [[45]4]. R. henanense subsp. lingbaoense
   (Rhl) is a rhododendron subspecies from Henan Province, endemic to
   Lingbao. It is distributed in the mountaintops or ridges above 2000 m
   elevation in Laoyacha [[46]5]. Rhl is a woody species with relatively
   weak natural reproduction ability. It is an ancient rhododendron tree
   that is a remnant of the natural forest in the Xiaoqinling Mountains
   [[47]6] and has scientific research and natural heritage value as a
   germplasm resource bank [[48]7].

   Global temperatures have increased by approximately 1.1℃ compared to
   pre-industrial times [[49]8]. Global warming has also increased the
   frequency of extremely hot days, seriously affecting the growth,
   survival and distribution of vegetation worldwide [[50]9]. Heat stress
   can seriously damage the normal physiological functions of
   rhododendrons, which can lead to the wilting of branches and leaves or
   the death of plants [[51]10]. It causes oxidative burst within cells,
   damaging organelles and cell membranes and leading to enzyme
   denaturation [[52]11]. To cope with temperature fluctuations, plants
   have developed various mechanisms such as activation of signal
   molecules (MAPK, CDPK, etc.), antioxidant enzymes [catalase (CAT),
   peroxidase (POD), superoxide dismutase (SOD), glutathione peroxidase,
   etc.], heat shock factors (TaHSFA6e, HSF3, etc.), heat shock proteins
   (HSPs; HSP26, HSP17, sHSPs, etc.), starch reduction, synthesis of amino
   acids and proteins, and expression of many stress-related genes (ABC
   transporter, pyrolyte disulfide isomerase, etc.) [[53]12, [54]13].
   Moreover, dehydration responsive element binding proteins 1 and 2 are
   important transcription factors (TFs) in the regulation of drought/heat
   stress tolerance [[55]14, [56]15]. Heat stress (HS) induces HSPs
   including HSP26, HSP17, and HSP70 and a complex transcription network
   including DREB2A and heat shock factors a1 (HSFA1s). Among them, HSFA1s
   is the master regulator of the HS response, which directly regulates
   important TFs including HSFA2, DREB2A, and HSFBs [[57]16]. Four main
   pathways play important roles in response to heat stress in plants;
   they include the heat stress factor–heat shock protein (HSF–HSP),
   calcium ion–calmodulin, hormone regulation, and reactive oxygen species
   (ROS) [[58]17–[59]19].

   Rhododendron plants exhibit certain heat resistance; however, most
   species can not withstand scorching sun and tolerate cold and humid
   environment. When the ambient temperature is higher than 35℃,
   rhododendron plants enter dormancy; strong light causes burning and
   yellowing of the leaves, and the whole plant may die [[60]20]. It was
   reported that with the increase of temperature, the plant height and
   diameter of Rhododendron. staghorn, Rhododendron. rivulare increased;
   the activity of SOD, contents of chlorophyll, malondialdehyde(MDA),
   soluble protein, and soluble sugar increased, whereas the relative
   water content and total biomass of leaves significantly decreased
   [[61]21]. At present, high temperature in summer has become the main
   factor limiting the widespread application of rhododendron in urban
   green spaces [[62]3].

   Currently, combined metabolome and transcriptome analyses are widely
   used in plant research and have become a powerful tool to study the
   genes participating in the biosynthesis of various secondary
   metabolites in higher plants [[63]22]. A previous study used
   metabolomic and transcriptomic analyses to assess the response of maize
   to high temperature stress (HTS) and revealed that the downregulation
   of genes encoding sucrose phosphate synthase, soluble starch synthase,
   and amyloclase was closely related to decreased starch content, whereas
   upregulation of genes encoding glutamate synthase and glutamate
   pyruvate transaminase led to the increase in protein content [[64]23].
   Transcriptome and metabolome analyses were performed on pitaya
   seedlings exposed to high temperature; 64 differentially expressed
   metabolites (DEMs), including sugars, organic acids, and amino acids,
   and 198 differentially expressed genes (DEGs) were screened. HuPR-1 was
   confirmed to play a significant role in the response of pitaya to HTS
   [[65]24]. Metabolome studies reported that in Arabidopsis thaliana,
   sucrose can replace proline as the main osmoprotectant of cell membrane
   under HTS, and the contents of proline and soluble sugar in metabolic
   network play a very important role in the response of plants to HTS
   [[66]25]. The physiological and biochemical responses of leaves from
   four rhododendron varieties under HTS were studied; soluble protein
   content exhibited a significant downward trend with the increase of
   temperature, and their physiological metabolism significantly changed
   [[67]26]. For example the contents of MDA, proline, reduced glutathione
   (GSH), superoxide dismutase (SOD), catalase (CAT), and peroxidase (POD)
   in Rhl leaves exhibited an overall trend of first increasing and then
   decreasing with the increasing temperature.In recent years, the studies
   on rhododendron are mainly focused on garden application [[68]27],
   tissue culture [[69]28], phytochemistry and pharmacology [[70]29],
   drought stress, and response to cold and HS [[71]30–[72]32]. The
   molecular mechanism of response of rhododendron to high temperature has
   not been thoroughly studied in China and abroad. As an endangered
   alpine plant, Rhl should be studied in detail; however, studies on Rhl
   are scarce. Recently, Rhl genome was assembled from scratch using
   PacBio long-read and Illumina short-read sequencing technologies
   [[73]33]. The C2H2 gene family of Rhl had been identified, and its
   expression has been analyzed under different stress treatments
   [[74]34]. The transcriptome and metabolome of Rhl under HTS have not
   been analyzed. Therefore, in this study, Rhl seedlings were treated at
   40℃ (RLH), 32℃ (RLM), and 24℃ (RLC), and the changes in transcriptome
   and metabolome were compared. The dynamic changes of the gene
   expression and metabolite accumulation in Rhl leaves at different
   temperatures were analyzed using RNA-seq, and liquid
   chromatography–tandem mass spectrometry (LC-MS/MS). qRT-PCR analysis
   was performed for the validation of transcriptome sequencing results.
   The results of this study can provide abundant genomic information for
   future research and help to clarify the molecular mechanisms of heat
   stress response of Rhl.

Materials and methods

Plant materials

   In October 2023, an appropriate amount of Rhl seeds was collected from
   the Xiaoqinling National Nature Reserve in western Henan Province. Rhl
   seeds were sterilized with 2% sodium hypochlorite for 10–15 min and
   washed 5 times with sterile waterFurther, the seeds were germinated
   using double filter paper method in culture dish. The seedlings were
   cultured under a 16 h light/8 h darkness photoperiod condition in a
   light incubator with the humus soil in the native territory
   (tepperature: 24℃, light intensity: 10,000 lx, and relative humidity:
   75%.) until the third true leaf stage. They were divided into three
   groups: control group (RLC), medium temperature group (RLM), and high
   temperature group (RLH). The seedlings in these groups were placed at
   24℃, 32℃, and 40℃ in a Light incubator (Bluepard, Shanghai Yiheng
   Technology Co., LTD, CN), respectively, for 72 h. Further, tender
   leaves were collected, freezen in liquid nitrogen, and stored at -80℃.
   For metabolome and transcriptome analyses, six and three replicates per
   group were used, respectively.

Metabolome analysis

Pre-processing

   An appropriate amount of leaves was accurately weighed in a tube
   (100 mg). To this, 600 µL methanol (Fisher Scientific, Loughborough,
   UK). was added and vortexed for 30 s. To this, a steel ball was added.
   The tube was placed in a tissue grinder for 60 s at 55 Hz, followed by
   ultrasonication for 15 min. Further, the solution was centrifuged at
   13,400 rpm at 4 °C for 10 min, and the supernatant was filtered through
   a 0.22-µm membrane. The filtrate was added to the detection bottle for
   LC-MS.

   The LC analysis was performed on a Vanquish UHPLC System (Thermo Fisher
   Scientific, USA). Chromatography was conducted using an ACQUITY UPLC^®
   HSS T3 (2.1 × 100 mm, 1.8 μm; Waters, Milford, MA, USA). The column:
   T3, the reverse column, column temperature: 40℃. The flow rate and
   injection volume were set at 0.3 mL/min and 2 µL, respectively. For
   LC-ESI(+)-MS analysis, the mobile phases consisted of 0.1% formic acid
   in acetonitrile (v/v) (B2) and 0.1% formic acid in water (v/v) (A2).
   Gradient program was as follows: 0–1 min, 8% B2(92% A2); 1–8 min, 8–98%
   B2(92–2% A2); 8–10 min, 98% B2(2% A2); 10–10.1 min, 98–8% B2(2%–-92%
   A2); 10.1–12 min, 8% B2(92% A2). For LC-ESI(-)-MS analysis, the mobile
   phase was acetonitrile (B3) and 5 mM ammonium formate (A3). The
   gradient elution procedure was the same as that for the cationic mode.
   MS detection of metabolites was performed on Orbitrap Exploris120
   (Thermo Fisher Scientific, USA) with ESI ion source. Simultaneous MS1
   and MS/MS (Full MS-ddMS2 mode, data-dependent MS/MS) acquisition was
   used.

Data analysis

   The original data were first converted into mzXML using MSConvert in
   ProteoWizard package software (v3.0.8789) [[75]35] and further
   processed using R XCMS (v3.12.0) for feature detection [[76]36]. The
   metabolites with relative standard deviation(RSD) > 30% in the quality
   control (QC) samples were filtrated. Subsequently, identification was
   conducted through precise MS data and matched with massbank [[77]37],
   HMDB [[78]38], and KEGG [[79]39].

   Unsupervised and supervised multivariate statistical analysis models
   [principal component analysis (PCA) and orthogonal partial least
   squares discriminant analysis (OPLS-DA), respectively] were used to
   differentiate the different groups using R ropls (v1.22.0) software
   [[80]40]. The steps of data analysis for 6 duplicate samples in each
   group (RLC, RLM, and RLH) were as follows: FastQC was used for raw data
   quality control; Trimmomatic for trimming and filtering; STAR was used
   to compare the data with the reference genome; featureCounts was used
   to quantify gene expression; DESeq2 was used for differential
   expression analysis, and the false discovery rate was controlled using
   the Benjamini-Hochberg program. To identify the DEMs, the screening
   criteria were fold change ≥ 2 or ≤ 0.5, P-value < 0.05, and variable
   importance in the project (VIP) ≥ 1 [[81]41]. The DEMs were subjected
   to KEGG pathway enrichment analysis using MetaboAnalyst
   ([82]www.metaboanalyst.ca) [[83]42].

Transcriptome sequencing analysis

   The RNA-Seq of the leaves from 9 samples (RLC, RLM, and RLH) was
   performed by Frasergen Biotechnology Co., Ltd. (Wuhan, China). After
   passing the three library tests, DNA nanoballs were prepared and
   mounted on sequencing chips. MGISEQ-T7(PE150) sequencing platform was
   used for sequencing. The library preparation and quality test were
   performed as follows. RNA sample was extracted, and RNA quality was
   evaluated using NanoDrop spectrophotometer and Agilent Bioanalyzer; The
   cDNA library was constructed using TruSeq RNA sample preparation kit;
   Qubit and Agilent high sensitivity DNA kits were used for quantitative
   analysis and quality detection of the library. Further, SOAPnuke
   software (v2.1.0) was used to filter raw reads to obtain high-quality
   clean reads. HISAT2 software was used to compare clean data with the
   Rhl reference genome [[84]33], and alignment quality was assessed using
   RSeQC (v3.0.1). For the known genes, QC, database annotation,
   quantitative expression, analysis of significant difference
   (P-value < 0.05), GO and KEGG functional enrichment analyses, mining of
   DEGs, and TF expression analysis were conducted. GO and KEGG enrichment
   analyses were conducted using Fisher’s exact test, and core genes
   within the GO/KEGG annotated gene set were identified using GSEA (Gene
   Set Enrichment Analysis, [85]http://gsea-msigdb.org/gsea/). For NR
   (NCBI non-redundant protein sequence), PFAM (protein family), and
   SWISSPROT (manually annotated and reviewed protein sequence database)
   annotations, a homology based approach was applied. BLAST was used for
   sequence comparison with the NR database.
   ([86]https://blast.ncbi.nlm.nih.gov/Blast.cgi). InterProScan was used
   to identify protein domains and functional sites, which aid in PFAM
   annotation([87]https://www.ebi.ac.uk/interpro/about/interproscan/).
   HMMER was used to search protein sequences in the PFAM database.
   ([88]http://hmmer.org/). UniProt was used to access SWISSPROT data and
   perform homology-based annotations([89]https://www.uniprot.org/). TFs
   were predicted using PlantTFDB 4.0. DESeq2 was used for the analysis of
   DEGs from each group of replicated samples. The analysis included
   standardization of counting data and Estimated dispersion. The DEGs
   were identified based on adjusted thresholds of P < 0.05 and log[2]
   fold change > 1. These steps ensured the robustness and reliability of
   the results, providing a clear understanding of the differences in gene
   expression between the groups.

Combined transcriptome and metabolome analyses

   Through KEGG pathway enrichment analysis([90]https://www.kegg.jp), the
   DEGs and DEMs were annotated to multiple metabolic pathways, and
   relevant metabolic pathways were screened and analyzed. In addition, R
   software was used to calculate the Pearson’s correlation coefficient of
   DEGs and DEMs and create the correlation heatmap. The DEGs and DEMs
   with P-value < 0.05 were selected to draw the network map using R
   version 3.5.1 software. The original transcriptome sequencing data have
   been submitted to the NCBI BioProject database with login number:
   PRJNA1141191.

Quantitative real-time PCR (qRT-PCR) analysis

   Total RNA was isolated using a Plant magnetic bead RNA extraction kit
   (Lingjun Biotechnology Inc., Shanghai, China). RNA concentration and
   purity were assessed using a NanoDrop-2000 (Thermo Fisher Scientific
   Inc., Waltham, USA). The first-strand cDNA was synthesized using Rever
   Tra™ kit (TOYOBO, Shiga, Japan), with 2 µg total RNA(DNA contamination
   had been removed), as per the manufacturer’s protocol. qRT-PCR analysis
   was performed using the CFX96™ real-time detection system (Bio-Rad
   Laboratories, Inc., Hercules, USA). All reactions were performed using
   the BestarTM qPCR MasterMix according to the manufacturer’s protocol
   (Cat. No. 2026789, DBI, DE), with three biological and three technical
   replicates for each sample. The primers of 9 genes and internal
   reference gene EF1α are given in Table [91]S1. The expression was
   calculated using 2^−△△CT method. SPSS 13.0 was used to analyze the
   results using ANOVA.

Results

Metabolome analysis

Sample quality control (QC) and data analysis

   PCA was performed using the datasets obtained from the RLC, RLM, and
   RLH samples The results indicated that the metabolites from three
   groups of samples treated at different temperatures were clearly
   separated in the score map, where the first principal component (PC1)
   and the second principal component (PC2) were separately plotted. For
   ESI^+ models, PC1 and PC2 accounted for 18.4% and 9.8% of the total
   change, respectively. For ESI^− models, PC1 and PC2 accounted for 17.2%
   and 9.8% of the total change, respectively (Fig. [92]1A and B). The
   positions of the three sample points on the PCA were relatively far,
   indicating significant differences among them. The QC samples were
   tightly clustered and located near the center of all samples,
   indicating good reproducibility of the experiment. The OPLS-DA plots
   were used to model the DEMs among RLC, RLM, and RLH. This model is
   stable and reliable and can be used to screen DEMs. At each
   temperature, both ESI modes could effectively separate the RLC, RLM,
   and RLH samples (Fig. [93]1C and D). These results indicated
   significant differences among the three treatments.

Fig. 1.

   [94]Fig. 1
   [95]Open in a new tab

   Principal component analysis (PCA) (QC sample) and orthogonal partial
   least squares discriminant analysis (OPLS-DA) scores of RLH group
   versus RLM group versus RLC group (A and C) Positive ion models. (B and
   D) Negative ion models

Analysis of differentially expressed metabolites (DEMs)

   A total of 78 DEMs were screened among the three treatments (Table
   [96]S2) with the criteria VIP value > 1 and P-value < 0.05. The DEMs
   were mainly distributed in fatty acyls group (12.8%), benzene and
   substituted derivatives (10.3%), carboxylic acids and derivatives
   (10.3%), and others (23.1%). A total of 76 DEMs were detected by
   comparing RLH vs. RLC (Fig. [97]2A). Among them, 40 and 36 DEMs
   exhibited increased and decreased levels, respectively; they included
   some amino acids namely L-proline, isoleucine, leucine, ornithine,
   methionine, and histidine (Fig. [98]2A, C). In RLM vs. RLC, the
   increase in amino acids content was insignificant. However, in RLH vs.
   RLC, it was significant by 2.11–5.33 times, indicating that HTS
   significantly increased the content of amino acids.

Fig. 2.

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

   The differentially expressed metabolites (DEMs) among three
   comparison.(A) The DEMs of up- (orange) and down-regulated (blue)
   metabolites in different comparison. (B) Venn diagram of DEMs. (C, D
   and E) Heatmap of DEMs. The complete linkage hierarchical cluster
   analysis (HCA) method was applied to reflect the normalized metabolite
   content in each sample, which is represented by a single row with
   ranging color chroma from red (high abundance) to blue (low abundance).
   C: RLC vs. RLH; D: RLH vs. RLM; E: RLC vs. RLM

   In total, 90 DEMs were detected in RLH vs. RLM. Among them, 33 and 57
   exhibited increased and decreased levels, respectively (Fig. [101]2A
   and D). Overall, 43 DEMs were detected in RLM vs. RLC. Among them, 29
   and 14 exhibited increased and decreased levels, respectively
   (Fig. [102]2A and E). In conclusion, the alterations in the metabolites
   in RLH vs. RLM and RLH vs. RLC were greater than those in RLM vs. RLC,
   and the metabolite content of Rhl under HTS(40℃) was significantly
   different compared with that under 35℃ and 24℃(Fig. [103]2C-E).

   RLH vs. RLC, RLH vs. RLM, and RLM vs. RLC had 27, 25, and 15 unique
   DEMs, respectively (Fig. [104]2B). Among the unique DEMs in RLM vs.
   RLC, the contents of L-rhamnono-1,4-lactone, (-)-bornesitol, and
   6-hydroxyhexanoic acid were significantly increased by 3.12–3.99 times
   (Tables [105]S3–[106]S5; P-value < 0.05). Among the unique DEMs in RLH
   vs. RLC, the contents of L-methionine, D-ornithine, isoquercetin,
   cyclic guanosine phosphate, and quercetin were significantly increased
   by 3.11–6.62 times (P-value < 0.05). Among the unique DEMs in RLH vs.
   RLM, the contents of 2-deoxystreptamine, 3-methyl-L-tyrosine, and
   bovinic acid were significantly increased by 2.21–2.64 times
   (P-value < 0.05).

KEGG pathway enrichment analysis

   In RLH vs. RLC, 76 DEMs were annotated to 50 pathways. Among them, the
   significantly enriched metabolic pathways were flavone and flavonol
   biosynthesis, D-amino acid and ABC transporters transport, and others
   (Fig. [107]3A). In RLH vs. RLM, 90 DEMs were enriched in 56 metabolic
   pathways. Among them, the significantly enriched metabolic pathways
   were linoleic acid metabolism, flavonoid and flavonol biosynthesis,
   tyrosine metabolism, ABC transport, and other pathways (Fig. [108]3B).
   In RLM vs. RLC, 43 DEMs were annotated to 29 metabolic pathways
   (Fig. [109]3C). Among the above pathways, the DEMs enriched in D-amino
   acid metabolism and ABC transporter pathway were more in number, and
   their P-values were relatively small. The impact value of flavonoid and
   flavonol biosynthesis pathway was larger (Fig. [110]3).

Fig. 3.

   [111]Fig. 3
   [112]Open in a new tab

   KEGG pathway enrichment analysis of DEMs in R. henanense subsp.
   lingbaoense after various temperature treatments (A) RLH vs. RLC. (B)
   RLH vs. RLM. (C) RLM vs. RLC. The dot color represents the P-value; the
   redder it is, the more significant the enrichment. The size of the dots
   represents the number of differential metabolites enriched in the
   pathway

   Among the DEMs with increased expression in RLH vs. RLC, two (quercetin
   and isoquercetin) were involved in flavonoid and flavonol biosynthesis
   pathway, and isoquercetin was the unique DEM (Fig. [113]S1A). Six DEMs
   (ornithine, L-leucine, L-histidine, L-proline, mannitol, and
   L-isoleucine) were involved in the D-amino acid metabolic pathway,
   among which L-methionine, D-ornithine, and cis-4-hydroxy-d-proline were
   the unique DEMs (Fig. [114]S1B). The results indicated that D-amino
   acid metabolism and flavonoid and flavonol biosynthesis pathway are
   affected under HTS, resulting in the upregulation of metabolites. The
   DEMs enriched in ABC transporter pathway mainly included phosphate and
   amino acid transporters, exhibited increased expression (Fig.
   [115]S1C). ABC transporter pathway is an important pathway for cell
   transport and may be an important pathway in the response of Rhl
   seedlings to HTS.

Transcriptome analysis

Quality evaluation

   A total of 66.798 Gb clean data (after QC) was obtained, and the
   average amount of clean data of each sample was 7.422 Gb. The
   percentage of Q30 bases was > 95.38%, and the GC content was
   45.76–46.16% (Table [116]S6). The results of 18,110, 10,157, 26,314,
   20,932, 358, and 18,480 genes have been successfully noted on the GO,
   KEGG, NR, PFAM, STRING, and SWISSPROT databases, respectively.

Detection of differentially expressed genes (DEGs)

   In total, 8450 DEGs were identified by comparison, of which 6110 DEGs
   were detected in RLH vs. RLC (3770 upregulated and 2340 downregulated);
   4429 DEGs were identified in RLH vs. RLM (2036 up-regulated and 2393
   down-regulated); and 6127 DEGs (3835 upregulated and 2292
   downregulated) were detected in RLM vs. RLC (Fig. [117]4A). Venn
   diagrams were drawn for the DEGs in the three comparison groups, with
   944 co-expressed genes (Fig. [118]4B). These results indicated that the
   expression levels of these DEGs may change during HTS.

Fig. 4.

   [119]Fig. 4
   [120]Open in a new tab

   The differentially expressed genes (DEGs) among three comparisonin. (A)
   The DEGs of up- (red) and down-regulated (blue) genes in different
   comparison. (B) Venn diagram of DEGs.3.2.3. Analysis of DEGs

   KEGG pathway enrichment analysis was performed on the DEGs
   (Fig. [121]5A-C). In RLH vs. RLC, DEGs were majorly enriched in
   photosynthesis, secondary metabolite biosynthesis, carbon metabolism,
   flavonoid biosynthesis, metabolic pathways, and other pathways.
   Significantly enriched pathways for the DEGs in RLH vs. RLM included
   metabolic pathways (P-value < 0.001), secondary metabolic biosynthesis
   (P-value < 0.001), plant hormone signal transduction (P-value < 0.001),
   etc. The enriched pathways for the DEGs in RLM vs. RLC mainly included
   metabolic pathway, signal pathway, carbon metabolism, secondary
   metabolic biosynthesis, plant hormone signal transduction, flavonoid
   biosynthesis, etc. These pathways in Rhl may be important in resisting
   HTS. GO pathway enrichment analysis revealed that 4095, 684, and 1724
   genes were significantly enriched in biological process (BP), cell
   composition (CC), and molecular function (MF), respectively (Fig.
   [122]S2). The most genes (1753) were enriched in Membrane
   (Fig. [123]5D). Moreover, 1346, 1334, 1249, and 1195 genes were mainly
   enriched in stimulus, membrane part, membrane intrinsic component, and
   membrane integral component, respectively.

Fig. 5.

   [124]Fig. 5
   [125]Open in a new tab

   KEGG and GO pathway enrichment analyses of DEGs in Rhl. KEGG pathway
   enrichment analysis of DEGs in (A) RLH vs. RLC. (B) RLH vs. RLM. (C)
   RLM vs. RLC. GO pathway enrichment analysis of DEGs in (D) RLH vs. RLC.
   The dot color represents the P-value; the redder it is, the more
   significant the enrichment. The size of the dots represents the number
   of differential genes enriched in the pathway

Mining of DEGs

   GO pathway enrichment analysis of DEGs in RLH vs. RLC revealed that
   several genes were enriched in the response to stimulus pathway in
   Biological process(BP), from which genes related to oxidative stress
   and HS were screened (Fig. [126]S2). KEGG pathway enrichment analysis
   revealed that photosynthesis and secondary metabolite biosynthesis were
   the significantly enriched pathways. Glycolysis/gluconeogenesis pathway
   can provide glucose for photosynthesis for secondary photochemical
   reaction. Analysis of the genes associated with these three pathways
   can help to reveal the mechanism of response of Rhl to HTS
   (Fig. [127]5A).

DEGs related to glutathione-S-transferase (GST)

   Under HTS, plants produce excessive ROS, and antioxidant enzymes play
   an important role in removing harmful ROS. Glutathione-S-transferase
   (GST) is a class of antioxidant enzymes involved in the metabolic
   pathway of glutathione and can play a role in protecting organisms from
   high temperatures. Overall, 13 genes related to GST were screened from
   the DEGs (Table [128]S7). Under HTS, except Rhe005696 and Rhe017231,
   the other 11 genes were upregulated.

DEGs related to HS response and photosynthesis

   Under HS, HSPs are produced to protect plants. In RLH and RLC, 43 DEGs
   were related to HSP (Table [129]S8), of which 35 were significantly
   downregulated and 8 were significantly upregulated. These results
   suggested that Rhl may inhibit the expression of genes encoding HSPs
   for self-protection in high temperature environments. It might be due
   to the tolerance of the plant to high temperatures. Photosynthesis is
   more sensitive to temperature. Compared with RLC, all 31 significant
   DEGs related to photosynthesis in RLH were upregulated (Table [130]S9).

Genes related to secondary metabolites and Glycolysis/gluconeogenic
biosynthesis

   Various secondary metabolites produced by secondary metabolic pathways
   help plants in coping with adverse environments. The top five genes
   with the most significant upregulation (328 in total) were selected
   (Table [131]S10). The most significantly upregulated gene was Rhe001827
   encoding chalcone and stilbene synthases that affect the production of
   flavonoids, which not only have antioxidant properties but also protect
   plants from heat damage. Second, an antioxidant enzyme POD encoded by
   Rhe016769 is an oxidoreductase in plant tissues and organs; it can
   eliminate hydrogen peroxide accumulated in the body and protect plants
   [[132]13]. Rhe005479 and Rhe012217 encode terpene synthase (TPS), and
   catalytic production of terpene compounds by them plays a significant
   role in plant stress tolerance [[133]43]. Glycolysis/ gluconeogenesis
   is closely related to photosynthesis. Overall, 40 DEGs related to
   glycolysis/gluconeogenesis were analyzed, most of which were
   significantly upregulated. The top five are shown in Table [134]S11.

Transcription factor (TF) expression analysis

   To investigate the expression of TFs in Rhl after HTS, stress-related
   TF family members such as bHLH, HSF, bZIP, ARF, NAC, MYB, AP2/ERF,
   WRKY, and C2H2 were analyzed. Overall, 3, 10, 17, 9, 2, 65, 50, 32, and
   12 genes of bHLH, HSF, bZIP, ARF, NAC, MYB, AP2/ERF, WRKY, and C2H2
   families were screened, respectively. The heatmap of these TF families
   (Fig. [135]6) revealed that compared with RLC, most members of bHLH,
   bZIP, ARF, NAC, MYB, AP2/ERF, and WRKY TF families were upregulated in
   RLM and RLH, which may be the reason for the higher number of
   upregulated genes in different comparison groups (RLH vs. RLC, RLH vs.
   RLM, and RLM vs. RLC). Moreover, it was observed that some members of
   the WRKY and AP2 families were upregulated under moderate heat stress
   (RLM), but downregulated under high heat stress (RLH). The TFs directly
   related to heat stress were Rhe007074 (HSFB4), Rhe012339 (HSFB3), and
   Rhe020232 (HSFAB). They were significantly upregulated in RLM vs. RLC
   and RLH vs. RLC. Rhe013859 (HSB2A) had the highest expression levels in
   RLM, and Rhe008336 (HSF30) and Rhe006303 (HFA4B) were significantly
   downregulated (See supplementary material1). These results demonstrated
   that complex transcriptional networks associated with heat shock
   response (HSR) are composed of many TF families.

Fig. 6.

   [136]Fig. 6
   [137]Open in a new tab

   Heatmap of differentially expressed Rhl transcription factors in the
   transcriptome. (A) Heatmap of 17 differentially expressed bZIP TFs. (B)
   10 differentially expressed HSF TFs. (C) 3 differentially expressed
   bHLH TFs. (D) 32 differentially expressed WRKY TFs. (E) 2
   differentially expressed NAC TFs. (F) 65 differentially expressed MYB
   TFs, (G) 50 differentially expressed AP2 TFs. (H) 12 differentially
   expressed C2H2 TFs. (I) 9 differentially expressed ARF TFs. Orange
   indicates high expression, and blue indicates low expression

Combined transcriptome and metabolome analyses

   Combined analysis of DEMs and DEGs among different samples was
   conducted. In RLH vs. RLC, RLH vs. RLM, and RLM vs. RLC, 45, 48, and 27
   DEM and DEG co-annotation pathways were observed (Supplement
   material[138]2-[139]1-[140]3). In three comparison groups, 13 distinct
   enrichment pathways of DEMs and DEGs were noted(Supplement
   material[141]2-[142]4). These results indicated that HTS may induce
   changes in the metabolites of these pathways. The flavonoid
   biosynthesis pathway and ABC transport pathway closely related to HTS
   were selected for in-depth analysis.

Flavonoid biosynthesis pathway in RLH vs. RLC

   Both transcriptome and metabolome analyses in this study indicated a
   significant upregulation of metabolites from flavonoid biosynthesis
   pathway. Therefore, the expression of genes and metabolites related to
   flavonoid biosynthesis under HTS was further investigated. Compared
   with RLC, the expression level of Rhe002465, Rhe006335, Rhe014544,
   Rhe011722, Rhe014879, Rhe016929, Rhe019664, Rhe019800, Rhe024744,
   Rhe025440, Rhe025485, Rhe026310, and Rhe001827 was significantly
   increased, whereas that of Rhe006334, Rhe022609, and Rhe002395 was
   significantly decreased in RLH (Fig. [143]7A). Compared with RLC,
   quercetin content was higher and kaempferol content was lowered in RLH
   (Fig. [144]7B). The correlation between DEGs and DEMs revealed that
   kaempferol was positively correlated with Rhe002395 and Rhe022609
   (correlation coefficients: 0.6761 and 0.8052, respectively,
   P-value < 0.001). Conversely, quercetin was negatively correlated with
   them (Fig. [145]7C and D).

Fig. 7.

   [146]Fig. 7
   [147]Open in a new tab

   Correlation analysis of Flavonoid biosynthesis pathway between
   transcriptome and metabolic group in Rhl high temperature (RLH) and
   control group(RLC). (A) Flavonoid biosynthesis pathway related DEGs
   analysis. (B) Differential metabolite expression pattern analysis. (C
   and D) Analysis of correlation between DEGs and DEMs. Note: (A, B, and
   C) *: P < 0.05; **: P < 0.01; ***: P < 0.001. (D) Correlation
   coefficient > 0.8. Orange or red indicates high expression, and blue
   indicates low expression

   p-Coumaroyl-CoA is catalyzed by CHS to form chalcone, followed by CHI
   to form naringen, F3H to form dihydrokaempferol, FLS to form
   kaempferol, F3H to form dihydroquercetin, and finally FLS to form
   quercetin (Fig. [148]8). High temperature induced significant
   upregulation of CHS, E5.5.1.6, and F3H, the three synthetases in the
   flavonoid biosynthesis pathway (the significantly upregulated
   regulatory genes of three enzymes are given in Fig. [149]8) and further
   induced significant accumulation of quercetin. In conclusion, high
   temperature decreased kaempferol but promoted quercetin synthesis,
   thereby enhancing the stress tolerance of Rhl under HTS.

Fig. 8.

   [150]Fig. 8
   [151]Open in a new tab

   Flavonoid biosynthesis pathway in R. henanense subsp. Lingbaoense. The
   heatmaps were drawn according to the FPKM values from the transcriptome
   datasets. Columns and rows in the heatmap represent samples and genes,
   respectively. The color scale indicates fold changes (Log2) in gene
   expression. chalcone synthase (CHS), flavanone 3-hydroxylase (F3H),
   flavonol synthase (FLS)

ABC transport metabolic pathways in RLH vs. RLC

   In RLH vs. RLC, 7 DEGs and 6 DEMs were identified in the ABC transport
   pathway. The 7 DEGs belong to the ABC transporter A3, ABCB1, and ABCC10
   subfamilies, respectively (Table [152]S12). Analysis of DEGs revealed
   that compared with RLC, the expression of Rhe007348, Rhe018087,
   Rhe021651, and Rhe021648 was significantly increased, whereas that of
   Rhe012338, Rhe023749, and Rhe029044 was significantly decreased in RLH
   (Fig. [153]9A, Table [154]S12). DEM abundance analysis revealed that
   compared with RLC, L-isoleucine, L-leucine, L-proline, and ornithine
   levels significantly increased (P-value < 0.05), and mannitol levels
   decreased in RLH (Fig. [155]9B). Rhe007348 was not highly correlated
   with other metabolites except mannitol. The remaining genes were highly
   correlated with DEMs in the ABC transport pathway. Rhe018087,
   Rhe021648, and Rhe021651 were negatively correlated with mannitol but
   positively correlated with other metabolites (Fig. [156]9C and D).

Fig. 9.

   [157]Fig. 9
   [158]Open in a new tab

   Correlation analysis transcriptome and metabolome data of ABC metabolic
   pathway in Rhl between group high temperature (RLH) and control (RLC)
   groups. (A) Relevant DEMs. (B) DEM expression. (C and D) Analysis of
   correlation between DEGs and DEMs. Note: (A, B, and C) *: P < 0.05; **:
   P < 0.01; ***: P < 0.001; (D) Correlation coefficient > 0.8. Red or
   orange color indicates high expression, and blue color indicates low
   expression

qRT-PCR analysis for the validation of transcriptome sequencing results

   To further verify the reliability of transcriptome data, nine genes
   (Rhe000027, Rhe002895, Rhe004212, Rhe007247|, Rhe012415, Rhe012485,
   Rhe013826, Rhe015229, and Rhe023454) were randomly selected from 944
   DEGs among the three comparison groups for qRT-PCR analysis, and the
   correlation coefficients were all greater than 0.84 (Table [159]S13).
   The results of qRT-PCR were consistent with the trend of transcriptome
   sequencing results, indicating the reliability of RNA-seq data.
   (Fig. [160]10).

Fig. 10.

   [161]Fig. 10
   [162]Open in a new tab

   Verification of transcriptome results using qRT-PCR. Note: Different
   letters indicate significant differences in different temperature
   groups (P-value < 0.05)

Discussion

Heat stress increases metabolism of amino acids and flavonoids in Rhl

   In the field of Chinese medicine, rhododendron is one of the research
   hotspots. It has rich chemical composition and mainly contains volatile
   oil, diterpene, triterpene, sesquiterpene, steroid, alkanes and their
   oxygen-containing derivatives, alcohols, flavonoids, and coumarins
   [[163]44, [164]45]. R. irroratum and R. decorum have medicinal value
   [[165]46, [166]47]. Flowers of R. decorum can be used as food in local
   folk customs [[167]48]. It was reported that total flavonoids in R.
   argentosa extract had a certain therapeutic effect in ischemic heart
   disease [[168]49]. In this study, 78 DEMs were detected in Rhl after
   HTS; they were enriched in fatty acyl groups, carboxylic acids and
   their derivatives, isoprene lipids, benzene and its substituted
   derivatives, organic oxygen compounds, terpene lactones, and phenols
   (Table [169]S2). Considering these metabolites in Rhl, it may have
   significant medicinal value.

   Plants respond to HTS through a series of intracellular metabolic
   reactions, such as starch production, amino acid and protein syntheses,
   and production of abscisic acid, antioxidants, and other protective
   molecules [[170]50, [171]51]. In this study, Rhl accumulated several
   amino acids and their derivatives including proline, histidine, and
   leucine, which can be used as osmotic protectants and molecular
   chaperones. Moreover, they play an important role in eliminating ROS,
   stabilizing the protein structure, and maintaining intracellular
   homeostasis [[172]52] (Table. [173]2). Therefore, accumulation of amino
   acids and upregulation of related metabolic pathways are important ways
   for plants to cope with stress. KEGG pathway enrichment analysis
   revealed that the D-amino acid metabolism pathway promoted the
   expression of methionine, ornithine, proline, etc., whereas the ABC
   transporter pathway increased the transport of phosphate and amino
   acids. This demonstrated that amino acid metabolism and expression of
   corresponding amino acids were important pathways for Rhl to respond to
   high temperature. This was consistent with a study reporting that
   Spinacia oleracea promoted amino acid metabolism under HTS [[174]53].
   Another important enriched KEGG pathway was flavonoid and flavonol
   biosynthesis. Flavonoid compounds are rich in hydroxyl groups and have
   antioxidant effects; for example, phenolic hydroxyl groups in quercetin
   and isoquercetin provide hydrogen atoms to react with free radicals
   [[175]54]. In this study, the unique flavonoid DEMs in RLH vs. RLC were
   the most abundant, among which quercetin and isoquercetin were
   upregulated after HTS. This indicated that flavonoids play an important
   role in the tolerance of Rhl to HTS.

Expression analysis of key genes and transcription factors related to heat
stress in Rhl

   In this study, after RNA-seq of Rhl samples, 94,451 genes were obtained
   and annotated, which were close to the previously reported data of R.
   hainanense, R. molle, and R. irroratum [[176]55, [177]56]. Overall,
   8450 DEGs were identified in three comparison groups. 9 genes were
   randomly selected to velidate the reliability of RNA-seq results using
   qRT-PCR (Fig. [178]2) [[179]57]. The results suggested that heat
   stimuli significantly affected the transcriptome profiles in Rhl and
   provided more meaningful information for further studies on heat
   tolerance in rhododendron. GO pathway enrichment analysis revealed that
   the DEGs were enriched in response to heat, membrane, response to
   stimulus, ROS, and abiotic stimuli, which were the top enriched GO
   terms in RLH vs. RLC (Fig. [180]5D). KEGG pathway enrichment analysis
   revealed that 1116 DEGs were significantly enriched in biosynthesis of
   secondary metabolites, flavonoid biosynthesis, plant hormone signal
   transduction, glycolysis/ gluconeogenesis, etc. (Fig. [181]5A–C). These
   results indicated that signal pathways related to HS were multiple and
   complex in Rhl, which is consistent with model plants such as
   Arabidopsis and rice [[182]58, [183]59]. These results were important
   for providing reference resources for molecular research in
   rhododendron plants in the future.

   High temperature can damage plant cell membrane system [[184]60]. GO
   pathway enrichment analysis revealed that the number of genes enriched
   in membrane was the highest. It could be predicted that the response of
   Rhl to HTS was through its influence on the cell membrane. HS can
   stimulate ROS production because it increases the levels of superoxide
   anion (O[2]^−), hydroxyl radical (HO·), and hydrogen peroxide
   (H[2]O[2]) in the mitochondria [[185]61, [186]62]. Excessive
   accumulation of these oxygen radicals results in cellular damage and
   death [[187]63]. However, plants have evolved mechanisms to balance the
   production and scavenging of ROS through enzymatic and nonenzymatic
   antioxidative processes [[188]63]. High antioxidant capacity, such as
   that obtained by enhancing the activities of SOD, POD, CAT, ascorbate
   peroxidase, and glutathione reductase, is conducive to improving the
   heat tolerance of plants [[189]64, [190]65]. In this study, 11 genes
   related to GST were screened for upregulation. Rhe016769 encoding POD
   was significantly upregulated. It is speculated that the high
   antioxidant capacity acquired by enhancing the activity of GST and POD
   is beneficial for eliminating ROS and improving heat tolerance in Rhl.

   The heat tolerance of rhododendron is mainly regulated by the
   physiological function of the photosynthetic system [[191]66].
   Photosynthesis is very sensitive to temperature. In this study, 31
   genes related to photosynthesis were significantly upregulated,
   indicating that high temperature improves photosynthesis in Rhl. Under
   stress conditions, the ability of plants to absorb and utilize light
   energy is greatly reduced, resulting in increased light inhibition
   [[192]67]. Under HTS, several protein subunits in the light and
   electron transport chains (such as Psb27 and Psb28) and coenzymes (such
   as CP43) mediate PSII repair [[193]68, [194]69]. In this study,
   photosynthesis-related genes were upregulated. The could be because
   photosynthesis in Rhl was improved under HTS, or photosynthesis-related
   genes were downregulated and photosynthesis was inhibited at first
   [[195]70]. Further, with the increased duration of stress, Psb28 and
   other mediated PSII repair, and the photosynthesis-related genes were
   upregulated. However, this should be further studied and verified.
   Furthermore, photosynthesis is closely related to
   glycolysis/gluconeogenesis, which can provide glucose for secondary
   photochemical reactions. This pathway can be used as a protective
   mechanism of cells under HTS [[196]71]. Most genes regulating
   glycolysis/gluconeogenesis-related proteins were significantly
   upregulated, and the regulation of these genes, photosynthesis related
   genes, antioxidant related genes, and HSP-related genes may alleviate
   HTS-induced damage to Rhl.

   HSP accumulation can enhance the heat tolerance of plants [[197]72].
   HSP90 plays a protective role in heat shock reaction [[198]73]. HSP70
   family proteins are expressed under stresses such as high temperature,
   cold, and stimulation by chemical reagents [[199]74]. Similarly, it was
   reported that with an increase in temperature, stress-associated genes
   including CfAPX1, CfAPX2, CfHSP11, CfHSP21, CfHSP70, CfHSFA1a, and
   CfHSFB4 were upregulated in Cryptomeria fortunei under HTS [[200]75].
   In this study, highly enriched HSPs were mostly downregulated under HTS
   (Table [201]S8) probably because of considerable difference between the
   high temperature to which the plants were exposed and optimal growth
   temperature. After the gene expression reached the maximum value, it
   exhibited a downward trend because of its intolerance to excessive
   temperature [[202]76]. It is always believed that HSPs can protect
   plants from damage. In this study, 35/43 DEGs related to HSP were
   significantly downregulated. Rhl may inhibit the expression of genes
   encoding HSPs for self-protection under high temperature environments,
   but further research is needed. It was concluded that flavonoids can
   help rice plants to resist high temperature and heat damage [[203]77].
   Significant upregulation of chalcone synthase encoded by Rhe001827
   would increase the yield of flavonoids and enhance the tolerance of Rhl
   to high temperature. This suggested that flavonoid biosynthesis is also
   an important pathway in Rhl in response to HTS, which is consistent
   with the metabolome analysis. This enhanced the understanding of
   molecular mechanism of response of Rhl to high temperature.

   TFs directly regulate the expression levels of stress-responsive genes
   and play a crucial role in the responses of plants to biotic and
   abiotic stresses [[204]78]. In this study, 9 TF families associated
   with stress were analyzed. The R2R3-MYB TF of arabidopsis AtMYB4 plays
   a vital role in regulating phenylpropanoid metabolism and also
   participates in abiotic stress responses [[205]79]. In this study,
   compared with RLC, 46 MYB family genes were upregulated (19
   downregulated) in RLH and RLM. AP2/ERF TF family is one of the major
   cold-stress-related TF families in plants. Together with other TF
   families including WRKY, bHLH, bZIP, and MYB, it enhances cold stress
   tolerance in plants [[206]80]. Compared with HLM, 46 AP2/ERF family TFs
   were upregulated and 4 were downregulated in this study. WRKY TFs are
   reported to be involved in response to various biotic or abiotic
   treatments [[207]81]. In A. thaliana, it is reported that WRKY71 is
   involved in regulating the ROS, SA, and jasmonic acid signaling
   pathways, which may promote the adaptation of plants to various abiotic
   stresses [[208]82]. In this study, compared with RLC, 31 WRKY family
   genes were upregulated (1 was downregulated) in RLH and RLM. Among
   them, upregulation of Rhe029954 (homologous to WRKY71) during HTS
   indicated its potential regulatory mechanism in HSR (Fig. [209]6).
   Moreover, it was observed that some members of the WRKY and AP2
   families were upregulated under moderate heat stress (RLM), but
   downregulated under high heat stress (RLH). This phenomenon may be
   related to the response mechanism of plants under different heat stress
   levels. This discrepancy might suggest that the transcription factors
   are responsive to heat stress in a dose-dependent manner, with some
   members of these families being more strongly upregulated at moderate
   stress levels and potentially reaching a saturation point or undergoing
   downregulation when the heat stress becomes more severe. This
   phenomenon could be indicative of a feedback mechanism or a cellular
   adjustment to prevent overexpression of stress-related genes under
   prolonged high stress conditions [[210]83–[211]84]. NAC TFs is a large
   family involved in various processes of plants, such as growth,
   development, and tolerance to stress [[212]85]. In this study, 2 NAC
   family genes were upregulated. HSFs play an important role in the
   signal transduction cascade, activating the genes of heat-induced
   response and acting as key TFs in the regulatory network of
   heat-induced response [[213]16]. A total of three HSF genes were
   identified in Rhl including Rhe007074, Rhe012339, and Rhe020232, which
   were significantly upregulated under HS. Therefore, these TF families
   might be the key regulators of HS response pathways in Rhl.

Heat stress increases metabolism and expression of genes of flavonoid
biosynthesis and ABC transport pathway in Rhl

   To obtain more comprehensive understanding of the metabolic processes
   and gene expression regulation mechanisms of Rhl under HTS, combined
   transcriptome and metabolome analysis was conducted. As a class of
   important plant secondary metabolites, flavonoid compounds play a
   crucial role in physiological and biochemical reactions, growth and
   development, stress, and stress responses in plants [[214]86]. In this
   study, KEGG pathway enrichment analysis revealed that flavonoid
   biosynthesis was significantly enriched after HTS (Fig. [215]5A), and 2
   DEMs and 19 DEGs related to flavonoid biosynthesis were obtained
   (Figs. [216]7 and [217]8).Among the 2 DEMs, kaempferol positively
   correlated with Rhe002395 and Rhe022609 (correlation coefficients:
   0.6761 and 0.8052, respectively; P-value < 0.001), but quercetin
   negatively correlated with them. Flavanone 3-hydroxylase (F3H) is one
   of the main factors in flavonoid biosynthesis. Heterologous expression
   of DoF3H in Escherichia coli led to higher cold and salt stress
   tolerance [[218]87]. Two F3H genes were reported in 19 DEGs, and HTS
   increased the expression levels of these two F3H genes (Fig. [219]8).
   These results indicated that these genes involved in the flavonoid
   synthesis pathway may enhance the high temperatures response of Rhl,
   and the increased synthesis of certain metabolites might improve
   adaptation to heat stress.

   The ABC transporters that have been extensively studied in plants are
   the ABCB, ABCC, and ABCG subgroups. In A. thaliana, the main function
   of the ABCB subfamily is to transport secondary metabolites and
   hormones and control stomata opening and closing. Moreover, ABCC
   subfamily is reported to be involved in glucose and chlorophyll
   metabolism in A. thaliana [[220]88, [221]89]. It was reported that
   CsABCC1, a major regulatory gene of the ABC transport pathway, may
   cause the accumulation of amino acids and carbohydrate metabolites and
   participate in the response to HTS [[222]90]. In this study, the
   correlation analysis of DEGs and DEMs under HS revealed that Rhe018087
   and Rhe021651 were highly positively correlated with L-isoleucine,
   L-leucine, and L-proline (Fig. [223]9). It can be concluded that the
   process of HTS causes changes in ABC transport metabolic pathway, and
   high temperature induces the increase in expression levels of Rhe018087
   and Rhe021651, and the accumulation of amino acids including L-proline.
   Rhe018087 and Rhe021651 may be involved in the response to HTS. It was
   observed that many high-temperature-responsive elements in Cajanus
   cajan, such as CcABCG7, CcABCG14, and CcABCG32 were significantly
   upregulated under HTS at 42 °C [[224]88]. ABCG transporters could
   positively affect plant’s response to adversity by transporting
   specific substances, however no ABCG-related DEG was identified in the
   ABC transport pathway of Rhl, which may be attributed to the
   differences in growth environment and species.

Conclusions

   In this study, 78 DEMs were screened by metabolome analysis after
   treating Rhl seedlings to various temperatures. The alterations in the
   expression of metabolites in RLH vs. RLM and RLH vs. RLC were greater
   than those in RLM vs. RLC. In total, 6.798 Gb clean data were acquired
   via RNA-seq, and 8450 DEGs were identified by comparison. KEGG pathway
   enrichment analysis revealed that DEGs in RLH vs. RLC were mainly
   enriched in photosynthesis, secondary metabolic biosynthesis, flavonoid
   biosynthesis, and other pathways. Mining of DEGs revealed that the
   expression of some genes related to heat stress function increased
   highly significantly, e.g., Gst-related gene Rhe008987 (Log2 2.3097),
   HSP70-related gene Rhe013774 (Log2 2.6519), CHS encoding gene Rhe001827
   (Log2 9.7935), POD encoding gene by Rhe016769 (Log2 8.1921), and
   ADH1-coding gene Rhe008676 (3.8776) associated with glycolysis/
   gluconeogenesis. This needs to be further studied. Correlation analysis
   of metabolome and transcriptome revealed that 12 metabolic pathways
   were enriched in three comparison groups. In the flavonoid biosynthesis
   pathway, high temperature inhibited the expression of Rhe022609 (FLS);
   decreased the expression of kaempferol; induced the expression of
   Rhe006335(FLS), Rhe016929(FL3H), and other genes; and accumulated
   quercetin, which improved the stress tolerance of Rhl under HTS. In the
   ABC transport pathway, high temperature induced the expression of
   Rhe018087 (ABCB10) and Rhe021651(ABCB19), resulting in the accumulation
   of L-proline and other metabolites. This suggested that Rhe018087 and
   Rhe021651 are involved in the response to HTS. This study laid a solid
   foundation for the detailed understanding of the physiological and
   molecular responses of Rhl to HTS.

Electronic supplementary material

   Below is the link to the electronic supplementary material.
   [225]Supplementary Material 1^ (24KB, xlsx)
   [226]Supplementary Material 2^ (69.5KB, xls)
   [227]Supplementary Material 3^ (73KB, xls)
   [228]Supplementary Material 4^ (56KB, xls)
   [229]Supplementary Material 5^ (48KB, xls)
   [230]Supplementary Material 6^ (54.4KB, docx)
   [231]Supplementary Material 7^ (590.8KB, docx)

Acknowledgements