Transcriptomic and Metabolomic Analysis of the Heat-Stress Response of Populus tomentosa Carr.

: Plants have evolved mechanisms of stress tolerance responses to heat stress. However, little is known about metabolic responses to heat stress in trees. In this study, we exposed Populus tomentosa Carr. to control (25 ◦ C) and heat stress (45 ◦ C) treatments and analyzed the metabolic and transcriptomic e ﬀ ects. Heat stress increased the cellular concentration of H 2 O 2 and the activities of antioxidant enzymes. The levels of proline, ra ﬃ nose, and melibiose were increased by heat stress, whereas those of pyruvate, fumarate, and myo-inositol were decreased. The expression levels of most genes ( PSB27 , PSB28 , LHCA5 , PETB , and PETC ) related to the light-harvesting complexes and photosynthetic electron transport system were downregulated by heat stress. Association analysis between key genes and altered metabolites indicated that glycolysis was enhanced, whereas the tricarboxylic acid (TCA) cycle was suppressed. The inositol phosphate; galactose; valine, leucine, and isoleucine; and arginine and proline metabolic pathways were signiﬁcantly a ﬀ ected by heat stress. In addition, several transcription factors, including HSFA2 , HSFA3 , HSFA9 , HSF4 , MYB27 , MYB4R1 , and bZIP60 were upregulated, whereas WRKY13 and WRKY50 were downregulated by heat stress. Interestingly, under heat stress, the expression of DREB1 , DREB2 , DREB2E , and DREB5 was dramatically upregulated at 12 h. Our results suggest that proline, ra ﬃ nose, melibiose, and several genes (e.g., PSB27 , LHCA5 , and PETB ) and transcription factors (e.g., HSFAs and DREBs) are involved in the response to heat stress in P. tomentosa .


Introduction
Climate change has increased the frequency of extreme high temperatures worldwide [1]. High temperature (HT) is a major environmental stress that limits plant growth, metabolism, and productivity. As sessile organisms, plants have developed complex and diverse systems to respond to high temperatures [2]. For example, plants use antioxidants to alleviate oxidative stress by scavenging excess reactive oxygen species (ROS) and change osmotic stress by accumulating osmoregulatory metabolites such as proline, glycine, betaine, and soluble sugars (e.g., sucrose, glucose, and galactose) [3]. In addition, plants can respond to heat stress by a wide range of transcriptional changes and produce heat-shock proteins (HSPs) to maintain cellular homeostasis [4] and improve heat tolerance. Therefore,

RNA Isolation and Illumina Sequencing
Total RNA was extracted from the CTRL and heat-stressed leaves (three biological replicates each) using the MiniBEST Plant RNA Extraction Kit (TaKaRa, Dalian, China). RNA concentration and integrity were evaluated using the Qubit RNA Assay Kit and Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA) and the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA), respectively. Sequencing libraries were generated using the NEBNext Ultra RNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA) in accordance with the manufacturer's instructions, and mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. cDNA synthesis was subsequently performed as described previously [21]. RNA sequencing was performed on the Illumina Hiseq 4000 platform. A Perl script was used to remove raw reads containing adapter sequences, poly-N, and low quality reads (>10% unknown nucleotides). The clean reads were mapped to the Populus trichocarpa Torr. & Gray reference genome ( [22]) using TopHat v2.0.12. HTSeq v0.6.1 was used to determine the read count of each gene, and gene expression levels were estimated as fragments per transcript kilobase per million fragments mapped (FPKM) values [23]. Differentially expressed genes (DEGs) were defined as those with an adjusted q-value <0.05 and a |log2 fold-change| >1, as determined by DESeq.

RNA Isolation and Illumina Sequencing
Total RNA was extracted from the CTRL and heat-stressed leaves (three biological replicates each) using the MiniBEST Plant RNA Extraction Kit (TaKaRa, Dalian, China). RNA concentration and integrity were evaluated using the Qubit RNA Assay Kit and Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA) and the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA), respectively. Sequencing libraries were generated using the NEBNext Ultra RNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA) in accordance with the manufacturer's instructions, and mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. cDNA synthesis was subsequently performed as described previously [21]. RNA sequencing was performed on the Illumina Hiseq 4000 platform. A Perl script was used to remove raw reads containing adapter sequences, poly-N, and low quality reads (>10% unknown nucleotides). The clean reads were mapped to the Populus trichocarpa Torr. & Gray reference genome ( [22]) using TopHat v2.0.12. HTSeq v0.6.1 was used to determine the read count of each gene, and gene expression levels were estimated as fragments per transcript kilobase per million fragments mapped (FPKM) values [23]. Differentially expressed genes (DEGs) were defined as those with an adjusted q-value < 0.05 and a |log2 fold-change| > 1, as determined by DESeq.
The sequencing data have been deposited in the Gene Expression Omnibus (GEO) database under accession number GSE124805.
Gene ontology (GO) enrichment analysis of DEGs was carried out with the GOseq R package [24]. All DEGs were assigned to pathways of the KEGG database, and their enrichment in DEGs was assessed using KOBAS software [25]. Corrected p-values (q-values) <0.05 were considered indicative of terms or pathways significantly enriched in DEGs.

Metabolite Preparation
Leaf tissue samples (0.06 g) were added to 2 mL EP tubes, soaked in 0.48 mL of extraction solvent (methanol:H 2 O, 3:1) with 24 µL (2 mg/mL) of adonitol as an internal standard. Extractions were uniformed in a ball mill for 4 min at 50 Hz, followed by ultrasonication for 5 min in ice water and centrifugation for 15 min at 13,000 rpm and 4 • C. Next, 0.35 mL of supernatant was transferred to a fresh 2 mL GC/MS glass flask and dried in a vacuum concentrator without heating.
The dried metabolites were added to 80 µL of methoxy amination hydrochloride (20 mg/mL in pyridine) and incubated for 30 min at 80 • C. Next, 100 µL of the BSTFA regent (1% TMCS, v/v) were added to the samples, followed by incubation for 1.5 h at 70 • C and blending for GC/TOF-MS analysis.

Metabolomic Analysis
GC/TOF-MS analysis was performed on an Agilent 7890 gas chromatograph, connected to a Pegasus HT time-of-flight mass spectrometer. The GC was equipped with a DB-5MS capillary column (5% diphenyl, 95% dimethyl polysiloxane). Helium was used as the carrier gas at a constant flow rate of 1 mL min −1 with a splitless injector (1 µL aliquot of the analyte). The temperature was 50 • C for the first minute, followed by an increase to 310 • C at 10 • C/min over 26 min, and a hold at 310 • C for 8 min. The temperatures of the injection, transfer line, and ion source were set at 280, 270, and 220 • C, respectively. The energy was −70 eV in electron-impact mode. MS data were acquired in full-scan mode at 50-500 m/z at a rate of 20 spectra per second after a solvent delay of 460 s. Chroma TOF 4.3X software and the LECO-Fiehn Rtx5 database (LECO Corporation, St. Joseph, MI, USA) were used for raw peak extraction, data baseline filtering, baseline calibration, peak alignment, deconvolution analysis, peak identification, and peak-area integration [26]. The retention time index (RI) method was used for peak identification, with an RI tolerance of 5000.
The variable importance in the projection (VIP) values of the first principal component >1 combined with Student's t-test (p < 0.05) were used to screen metabolites (each treatment involved at least six biological replicates).

Scanning Electron Microscopy
Small strips (approximately 0.5 × 0.5 cm) were sheared from the areas between the margin and midrib of fresh leaves, rapidly fixed in 2.5% glutaraldehyde at 4 • C for 24 h, and dehydrated in an ethanol series (30% to 100%). The samples were subjected to automatic critical point drying (Leica, Wetzlar, Germany). The dried samples were mounted in metal stubs and sputter coated (Bal-Tec, Balzers, Switzerland) with a thin conductive film of gold. The coated samples were examined and photographed using scanning electron microscopy (SEM) (Hitachi High-Technologies Corporation, Tokyo, Japan) at 15.0 kV [27]. The average area of stomatal opening was calculated from at least 40 stomata using ImageJ software (V1.51; National Institutes of Health, Bethesda, MD, USA).

Assay of Antioxidant Enzyme Activity
Fresh leaves were sampled and ground to a fine powder in liquid nitrogen. Samples (0.1 g) were mixed with distilled water (1 mL) in an ice bath, centrifuged at 8000× g for 10 min at 4 • C, and the supernatants were used for further experiments. The activities of superoxide dismutase (SOD), catalase (CAT), and peroxidase (POD) were measured using commercial kits (Suzhou Comin Bioengineering, Suzhou, China). Each sample comprised three to five leaves to ensure a sufficient quantity of leaf tissue, and three biological replicates per treatment were performed.
SOD activity was measured using the nitro-blue tetrazolium (NBT) method [28]. SOD can scavenge the superoxide anion radical (O 2− ), which is generated from xanthine morpholine by xanthine oxidase, and O 2− reduces NBT, as indicated by observing absorbance at 560 nm after 30 min. CAT and POD activity assays were performed in a water bath at 25 • C for 10 min and the absorbance of H 2 O 2 was read at 240 and 470 nm, respectively.

Soluble Sugar and Soluble Protein Assays
Samples (0.1 g) were mixed with 1 mL distilled water in plastic centrifuge tubes and transferred to a water bath for 10 min at 95 • C. After cooling, the samples were centrifuged at 8000× g at 25 • C for 10 min. The supernatants were transferred to 10 mL tubes and the volume equalized to 10 mL with distilled water. The soluble sugar concentration was determined by anthrone colorimetry using a commercial kit (Suzhou Comin Bioengineering, Suzhou, China).
Samples (0.1 g) were mixed with 1 mL distilled water in an ice bath, and centrifuged at 8000× g and 4 • C for 10 min. The soluble protein concentration was determined using the Coomassie Brilliant Blue method with a commercial kit (Suzhou Comin Bioengineering). Each treatment involved at least three biological replicates.

H 2 O 2 Assay
Samples (0.1 g) were mixed with 50 mmol·L −1 phosphate buffer (0.9 mL) in an ice bath, and centrifuged at 10,000× g for 10 min at 4 • C. The H 2 O 2 concentration was measured using a colorimetric method via a commercial kit (Nanjing Jiancheng Bioengineering Institute, Nanjing, China). Each treatment involved at least three biological replicates.

Real-Time Reverse Transcription PCR
To evaluate gene expression levels, we carried out quantitative reverse transcription-polymerase chain reaction (qRT-PCR). RNA samples (containing approximately 1 µg total RNA) were treated with gDNA Eraser to eliminate any contaminant gDNA, and RNA solution (20 µL) was subjected to reverse transcription with HiScript II Q Select RT SuperMix for qPCR (+gDNA wiper) following the manufacturer's instructions (Vazyme, Nanjing, China). Primer 5.0 was used to design gene-specific primers, and actin2 was used as the internal reference gene. qRT-PCR was performed using a Bio-Rad CFX96 Real-Time System (Bio-Rad, Hercules, CA, USA) with the ChamQ SYBR qPCR Master Mix (Perfect Real Time) (Vazyme, Nanjing, China) commercial kit. The qRT-PCR conditions comprised an initial denaturation at 94 • C for 30 s, followed by 40 cycles of 95 • C for 10 s, 60 • C for 30 s, 95 • C for 15 s, 60 • C for 60 s, and a final 15 s at 95 • C. The 2 −∆∆Ct comparative threshold cycle (Ct) method was used to evaluate gene expression levels. Bio-Rad CFX Manager v1.6.541.1028 software (Bio-Rad, Hercules, CA, USA) was used to calculate Ct values, and each reaction involved three biological replicates.

Statistical Analysis
The antioxidant enzyme activity, soluble sugar, soluble protein, H 2 O 2 , and qRT-PCR experiments were all repeated three times in our study. In stomatal variation experiment, we calculated 40 stomata randomly for each group. The experimental data between control and treatments were analyzed based on Student's t-test at p less than 0.05 (SPSS 18.0 software for Windows) (SPSS, Chicago, IL, USA) [29].

Stomatal Variation and Physiological Changes under Heat Stress
We used SEM to observe stomatal closure and opening. There was a marked difference in stomatal opening between 25 • C and 45 • C (Figure 2a  Compared to the control, the total soluble sugar per unit leaf fresh weight was decreased by 32.11% in heat stress (p < 0.01) (Figure 3a), whereas the H2O2 concentration increased by about 16.22% (0.01 < p < 0.05) (Figure 3b). However, the soluble protein content did not differ significantly between the control and heat-stressed plants (Figure 3c). In addition, compared to the control, SOD and CAT activity were increased 130.51% (P < 0.01) and 43.52% (0.01 < P < 0.05), respectively, by heat stress  Compared to the control, the total soluble sugar per unit leaf fresh weight was decreased by 32.11% in heat stress (p < 0.01) (Figure 3a), whereas the H 2 O 2 concentration increased by about 16.22% (0.01 < p < 0.05) ( Figure 3b). However, the soluble protein content did not differ significantly between the control and heat-stressed plants ( Figure 3c). In addition, compared to the control, SOD and CAT activity were increased 130.51% (P < 0.01) and 43.52% (0.01 < P < 0.05), respectively, by heat stress (Figure 3d,e). Nonetheless, there was no significant difference in POD activity between the control and heat-stressed plants (Figure 3f).

Transcriptomic Changes in Response to Heat Stress
Using the Illumina HiSeq 4000 platform, we generated 73,405,830 and 73,497,212 raw reads from control and heat-stressed plants, respectively. After quality checks and data cleaning, we obtained 71,828,148 (control) and 70,418,198 (heat stress) clean reads (Table S1). The clean datasets were mapped to the P. trichocarpa genome; 87.60% (control) and 83.82% (heat stress) of the sequences could be aligned with the reference genome.
A total of 7,738 DEGs between the control and heat-stressed plants were found. Among them, 3,888 DEGs were upregulated and 3,850 DEGs were downregulated in the heat-stressed relative to the control plants ( Figure S1). The GO terms response to abiotic stimulus, response to heat, membrane-bounded organelle, chloroplast thylakoid, chloroplast stroma, and photosynthesis membrane ( Figure S2a, arrows) were significantly enriched. The pyruvate metabolism, photosynthesis, pentose phosphate pathway, glycolysis/gluconeogenesis, arginine and proline metabolism, carbon metabolism, and carbon fixation in photosynthetic organisms Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were significantly enriched in DEGs ( Figure  S2b, arrows). Some important DEGs and were listed in Table S2, and DEGs more than 1.5-log2 fold change were listed in Table S3.

Genes Associated with Photosynthesis and Respiration
Photosynthesis is a heat-sensitive process. Accordingly, we performed KEGG and GO enrichment analyses to identify genes associated with the photosynthetic electron transport system. Compared to the control, the majority of DEGs related to photosynthesis electron flow were significantly downregulated by heat stress (Figure 4a), including nine DEGs related to photosystem

Transcriptomic Changes in Response to Heat Stress
Using the Illumina HiSeq 4000 platform, we generated 73,405,830 and 73,497,212 raw reads from control and heat-stressed plants, respectively. After quality checks and data cleaning, we obtained 71,828,148 (control) and 70,418,198 (heat stress) clean reads (Table S1). The clean datasets were mapped to the P. trichocarpa genome; 87.60% (control) and 83.82% (heat stress) of the sequences could be aligned with the reference genome.
A total of 7738 DEGs between the control and heat-stressed plants were found. Among them, 3888 DEGs were upregulated and 3850 DEGs were downregulated in the heat-stressed relative to the control plants ( Figure S1). The GO terms response to abiotic stimulus, response to heat, membrane-bounded organelle, chloroplast thylakoid, chloroplast stroma, and photosynthesis membrane ( Figure S2a, arrows) were significantly enriched. The pyruvate metabolism, photosynthesis, pentose phosphate pathway, glycolysis/gluconeogenesis, arginine and proline metabolism, carbon metabolism, and carbon fixation in photosynthetic organisms Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were significantly enriched in DEGs ( Figure S2b, arrows). Some important DEGs and were listed in Table S2, and DEGs more than 1.5-log 2 fold change were listed in Table S3.

Genes Associated with Photosynthesis and Respiration
Photosynthesis is a heat-sensitive process. Accordingly, we performed KEGG and GO enrichment analyses to identify genes associated with the photosynthetic electron transport system. Compared to the control, the majority of DEGs related to photosynthesis electron flow were significantly downregulated by heat stress (Figure 4a), including nine DEGs related to photosystem I (PSI), seven related to photosystem II (PSII), two DEGs encoding cytochrome b6f complex, two DEGs related to light harvesting complex II (LHC II), and four DEGs related to light harvesting complex I (LHC I). Among them, the expression of PSB27, PQL3, and PPD3 (PS II); PSAD, PSAG1, PSAG2, and PSAK (PS I); LHCB5 (LHC II), LHCA3, LHCA4, LHCA5, and LHCA6 (LHC I); and PETC were all 2-log-fold lower under heat stress (Figure 4a). The expression of PSB28 and PETB was also reduced under heat stress. Therefore, PS II and PS I were suppressed by heat stress. GO analysis showed that four DEGs related to NADH dehydrogenase activity and PETH1 and PETH2 (related to ferredoxin-NADP+ reductase activity) were downregulated by heat stress (Figure 4a). Additionally, four genes related to ATP synthesis-coupled electron transport were downregulated, and one (POPTR_0004s20010) was upregulated, by heat stress. Thus, the photosynthetic light reaction was inhibited by heat stress.
Regarding oxidative phosphorylation, compared to the control, 11 genes related to NADH dehydrogenase were downregulated by heat stress. NADUAFB1, for example, was 2.85-log-fold downregulated. Two DEGs (COX5B and COX6B) encoding cytochrome c oxidase were downregulated, as were eight DEGs related to F-type ATPase, indicating suppression of oxidative phosphorylation (Figure 4b).

Integrated Analysis of Osmotic and Resistant Substance Metabolism-associated Gene Expression and Metabolite Levels
The following KEGG pathways were enriched in DEGs under heat stress: valine, leucine, and isoleucine biosynthesis; galactose metabolism; and glycine, serine, and threonine metabolism (Figure 5a).
In plants, carbohydrates and proline play key roles in the osmoregulatory response to stress. The levels of several carbohydrate conjugates (such as myo-inositol, galactinol, raffinose, and melibiose) and proline were significantly increased by heat stress (Figure 6). The genes encoding 1-L-myo-inositol-1-phosphate synthase and inositol-1-monophosphatase (in the inositol phosphate and galactose metabolism pathways) were significantly downregulated by heat stress (Figure 6). This was associated with reduced myo-inositol levels under heat stress (Figures 5b and 6). Two DEGs related to galactinol (a metabolite of myo-inositol) synthesis and five DEGs related to raffinose synthesis were upregulated by heat stress. However, the melibiose level increased despite the downregulation of two related DEGs ( Figure 6). The gene encoding 1-pyrroline-5-carboxylate synthetase (in the arginine and proline metabolism pathway) was upregulated and the proline content was increased 6.03-log-fold by heat stress (Figures 5b and 7). The genes encoding phosphofructokinase, phosphoglyceratekinase, enolase, and pyruvatekinase (involved in glycolysis/gluconeogenesis) were upregulated by heat stress (Figure 7); however, two genes encoding phosphoglucomutase were downregulated. The genes encoding pyruvate dehydrogenase E1 component subunit beta-1, pyruvate dehydrogenase E1 component subunit beta-3, dihydrolipoamide acetyltransferase, isocitrate dehydrogenase, succinate-CoA ligase, fumarate hydratase, and ATP citrate lyase (TCA cycle) were downregulated, and that encoding aconitate hydratase was upregulated, by heat stress. This is consistent with the decreased pyruvate and fumaric acid levels under heat stress (Figures 5b and 7). Finally, the valine, leucine, and isoleucine levels were significantly increased, and the gene encoding acetolactate synthase was upregulated, by heat stress (Figures 5b and 7).

Integrated Analysis of Osmotic and Resistant Substance Metabolism-associated Gene Expression and Metabolite Levels
The following KEGG pathways were enriched in DEGs under heat stress: valine, leucine, and isoleucine biosynthesis; galactose metabolism; and glycine, serine, and threonine metabolism ( Figure  5a). An association analysis between metabolite levels and gene expression detected seven integrative pathways ( Figure 6; Figure 7): inositol phosphate metabolism; galactose metabolism; arginine and proline biosynthesis metabolism; glycine, serine, and threonine metabolism; glycolysis/gluconeogenesis; citrate cycle (TCA); and valine, leucine, and isoleucine biosynthesis.   In plants, carbohydrates and proline play key roles in the osmoregulatory response to stress. The levels of several carbohydrate conjugates (such as myo-inositol, galactinol, raffinose, and melibiose) and proline were significantly increased by heat stress (Figure 6). The genes encoding 1-L-myoinositol-1-phosphate synthase and inositol-1-monophosphatase (in the inositol phosphate and galactose metabolism pathways) were significantly downregulated by heat stress (Figure 6). This was associated with reduced myo-inositol levels under heat stress (Figures 5b and 6). Two DEGs related to galactinol (a metabolite of myo-inositol) synthesis and five DEGs related to raffinose synthesis were upregulated by heat stress. However, the melibiose level increased despite the downregulation of two related DEGs ( Figure 6). The gene encoding 1-pyrroline-5-carboxylate synthetase (in the arginine and proline metabolism pathway) was upregulated and the proline content was increased 6.03-log-fold by heat stress (Figures 5b and 7). The genes encoding phosphofructokinase, phosphoglyceratekinase, enolase, and pyruvatekinase (involved in glycolysis/gluconeogenesis) were upregulated by heat stress (Figure 7); however, two genes encoding phosphoglucomutase were Figure 7. Associated pathway analysis of metabolite levels and enzyme activities. Arginine and proline biosynthesis metabolism; glycine, serine, and threonine metabolism; glycolysis/gluconeogenesis; citrate cycle (TCA); and valine, leucine, and isoleucine biosynthesis pathways. Metabolites whose levels did not differ were omitted. Colors at right indicate the log 2 expression ratio [log 2 (heat/control_ mean)], and at left, the FPKM value of the control. (Levels 1-8, FPKM values of 0-10, 10-20, 20-40, 40-80, 80-160, 160-320, 320-640, and 640+, respectively). Large icons, metabolites; red, upregulation; blue, downregulation.

Transcription Factors Involved in Heat Treatment
A total of 553 DEGs were identified as transcription factors (TFs) of 68 TF families (Table S4). We examined five of these TF families: HSFs, MYBs, WRKYs, BZIPs, and C2H2s.
The expression levels of HSFA9, HSFA3, HSFA2, and HSF4 were significantly upregulated by heat stress. In contrast, that of HSFB4 was significantly downregulated. Compared with the control, MYB27 and MYB4R1 were upregulated by heat stress. The genes encoding the WRKY50 and WRKY13 TFs were downregulated by heat stress. Notably, the gene encoding the C2H2 and C2HC zinc finger superfamily protein (STOP2) was upregulated by heat stress, as was bzip60 (Figure 8).
The expression levels of HSFA9, HSFA3, HSFA2, and HSF4 were significantly upregulated by heat stress. In contrast, that of HSFB4 was significantly downregulated. Compared with the control, MYB27 and MYB4R1 were upregulated by heat stress. The genes encoding the WRKY50 and WRKY13 TFs were downregulated by heat stress. Notably, the gene encoding the C2H2 and C2HC zinc finger superfamily protein (STOP2) was upregulated by heat stress, as was bzip60 (Figure 8). The responses of P. tomentosa to heat stress are summarized in Table 1. The responses of P. tomentosa to heat stress are summarized in Table 1. Table 1. Summary of responses in P. tomentosa subjected to heat stress.

Effects of Heat Stress on Antioxidant Activity
Heat stress induces the production of reactive oxygen species (ROS) in plants [30]. However, if ROS production passes a threshold value, it can lead to oxidative damage. Thus, the oxidantantioxidant balance is critical for plant survival. For example, corresponding, plants can increase their production of antioxidant enzymes such as SOD, POD, and CAT, which scavenge and remove ROS [31,32]. In this study, the production of H2O2 was increased under heat stress (45 °C), accompanied by a significant increase in the activity of SOD and CAT. Similar phenomena have been reported in other woody plants [33]. Nevertheless, a previous study in poplar reported no significant difference in the H2O2 level between ambient temperature (25 °C) and high temperature (30 °C), and SOD and CAT activity were decreased at 30 °C [18]. Because poplar is heat tolerant, the decreased activities of SOD and CAT suggest that exposure to 30 °C does not impose heat stress. However, exposure to 45 °C increased the activities of antioxidant enzymes, suggesting cellular damage and enhanced scavenging of ROS.

Effects of Heat Stress on Photosynthesis
Photosynthesis is suppressed at temperatures outside the optimal range [34]. Heat stress affects the light-harvesting systems, electron transport, NADPH and ATP synthesis, photosynthetic carbon

Effects of Heat Stress on Antioxidant Activity
Heat stress induces the production of reactive oxygen species (ROS) in plants [30]. However, if ROS production passes a threshold value, it can lead to oxidative damage. Thus, the oxidant-antioxidant balance is critical for plant survival. For example, corresponding, plants can increase their production of antioxidant enzymes such as SOD, POD, and CAT, which scavenge and remove ROS [31,32]. In this study, the production of H 2 O 2 was increased under heat stress (45 • C), accompanied by a significant increase in the activity of SOD and CAT. Similar phenomena have been reported in other woody plants [33]. Nevertheless, a previous study in poplar reported no significant difference in the H 2 O 2 level between ambient temperature (25 • C) and high temperature (30 • C), and SOD and CAT activity were decreased at 30 • C [18]. Because poplar is heat tolerant, the decreased activities of SOD and CAT suggest that exposure to 30 • C does not impose heat stress. However, exposure to 45 • C increased the activities of antioxidant enzymes, suggesting cellular damage and enhanced scavenging of ROS.

Effects of Heat Stress on Photosynthesis
Photosynthesis is suppressed at temperatures outside the optimal range [34]. Heat stress affects the light-harvesting systems, electron transport, NADPH and ATP synthesis, photosynthetic carbon cycle, and the utilization of assimilates [35,36]. Damage to any of these components can reduce photosynthetic capacity. PS II is the most heat-sensitive component of the photosynthetic apparatus [16]. Several protein subunits (e.g., Psb27 and Psb28) and cofactors (e.g., CP43) in photosynthetic electron transport mediate PSII repair under heat stress [37,38]. Light-harvesting-related proteins (LHCs) are important for light harvesting and photoprotection, and the Lhca1-6 and Lhcb1-8 LHC sub-classes are associated with photoprotection of PSI and PSII, respectively [39]. In this study, several genes related to the photosynthetic electron transport system and ATP synthase were downregulated under heat stress. These results indicate that photosynthetic electron transport and photophosphorylation are suppressed by heat stress. In particular, PSB27, LHCA3, LHCA5, and LHCA6 were downregulated eight-fold under heat stress, suggesting inhibition of photoprotection of PSII and PSI. Cytochrome b6f also mediates PSII repair by driving ATP synthesis in chloroplasts [16], and PETB and PETC (encoding the cytochrome b6f complex) were significantly downregulated by heat stress, suggesting a negative effect on PSII repair. Overall, these results suggest that PSII and PSI were damaged at 45 • C, resulting in suppression of the photosynthetic machinery in P. tomentosa.

Effects of Heat Stress on Respiration
Plant respiration generally increases with temperature [40]. Stomatal regulation can increase transpiration to cool the leaf in response to elevated temperatures; in contrast, excessive heat induces stomatal closure, resulting in reduced photosynthesis [41]. In this study, the size of the stomatal opening decreased significantly at 45 • C, suggesting a negative effect on gas exchange in leaves. Respiration comprises glycolysis, the TCA cycle, and the mitochondrial electron transport chain (ETC). These pathways can produce energy equivalents and carbon skeletons for the biosynthesis of metabolites [42]. In glycolysis and TCA cycle, high temperature can increase the levels of sucrose and pyruvate but decrease those of sugar phosphates and some organic acids [43]. The TCA cycle also supports ATP synthesis, and is embedded in a wider metabolic network [44]. Generally, the levels of TCA cycle intermediates, such as pyruvic acid, fumarate, malate, and citrate, are affected by heat stress [45]. The ETC transfers electrons from NADH to oxygen (the ultimate electron acceptor), and NADH ubiquinone oxidoreductase in the mitochondrial ETC is thermolabile [40].
In this study, most of the glycolysis pathway-related DEGs were upregulated by heat stress, suggesting enhanced consumption of glucose and reduction of carbohydrate. These genes encoding ATP citrate lyase, isocitrate dehydrogenase, and fumarate hydratase (involved in the TCA cycle) were downregulated by heat stress. Moreover, the levels of TCA cycle metabolites, such as pyruvate and fumarate, were reduced by heat stress. These results suggest that heat suppressed the TCA cycle. Some genes related to respiratory electron transfer and oxidative phosphorylation, including NADUFA1, NADUFA5, NADUFB1, COX5B, and COX6B, were significantly downregulated by heat stress. In addition, all ATP citrate lyases were downregulated. Therefore, heat stress decreased respiratory electron transfer and oxidative phosphorylation, resulting in reduced oxidative respiration.

Osmotic and Resistant Substance Metabolism in Response to Heat Stress
The plant heat-stress response involves adjustment of cellular osmosis via biosynthesis of osmoprotectants. Galactinol-and-raffinose synthase is a key enzyme in the synthesis of raffinose-family oligosaccharides, which function as osmoprotectants in plant cells [46,47]. The levels of raffinose and melibiose were increased by heat stress in this study, and the genes encoding galactinol and raffinose synthase were significantly upregulated, suggesting their importance in the response to heat stress.
Proline accumulation is a common physiological response to heat stresses, and has negative effects on plants during heat stress [48]. However, proline can also function as a metabolic signal of activation of anti-stress response system, including ROS detoxification, induction of the specific gene expression, and protection of membrane integrity [49]. In this study, the proline level was increased by heat stress, and the important proline synthesis gene P5CS (POPTR_0006s04100) was upregulated. Therefore, P. tomentosa responds to heat stress by upregulating P5CS, resulting in increased proline levels.
The above results indicate that exposure to 45 • C imposed osmotic stress on P. tomentosa, which affected sugar biosynthesis, remobilization, or redistribution, and caused increased accumulation of osmolytes, resulting in osmotic adaptation. In addition, the comprehensive changes in sugars and organic acid might also contribute to heat stress defenses.

Transcription Factors Involved in the Response to Heat Stress
HSPs function as molecular chaperones to prevent accumulation of denatured or misfolded proteins and facilitate protein reactivation after stresses [2]. These HSPs include HSP100, HSP90, HSP70, HSP60, and small HSPs [2], of which HSP70 is the most abundant. The HSP90s are involved in the maintenance of phenotypic plasticity and developmental stability [50]. Hsp70 and Hsp90-6 were upregulated 7549-and 8.5-fold, respectively, after 5 h at 45 • C. qRT-PCR showed that Hsp70 and Hsp90-6 were upregulated after 2, 5, and 12 h at 45 • C. Thus, Hsp70 and Hsp90-6 are involved in the maintenance of heat-stress tolerance in P. tomentosa.
Complex transcriptional networks (e.g., involving HSFs, DREBs, WRKYs, and bZIPs) are involved in the response to heat stress [51]. HSFs regulate the expression of HSPs in response to heat stress; e.g., HsfA2 induces the expression of HSP70 and HSP101 [52] and protects against oxidative damage to organelles under heat stress [53]. HsfA3 confers increased thermotolerance and its expression is directly regulated by DREB2s under heat stress [54]. In poplar, DREB2, DREB2C, and HsfA3 were significantly upregulated under heat stress [16]. Here, HsfA2, HsfA3, HsfA4, HsfA9, DREB2, and DREB2E were significantly upregulated by heat stress. We speculate that the HSFAs and DREB2s TFs induce the expression of HSPs and are integrated into the complex stress-response networks that mediate stress tolerance in P. tomentosa.
Other TFs, such as bZIPs, MYBs, and WRKYs, are involved in stress responses. bZIPs are endoplasmic reticulum (ER) stress sensors in plants that regulate abscisic acid (ABA) and stress signaling, and contribute to stress tolerance [55]. Members of the MYB and WRKY families participate in plant responses to heat stress. For example, WRKY25 increases the expression of HsfA2, HsfB1, HsfB2a, and Hsp101, and WRKY39 enhances thermotolerance [56,57]. In this study, we identified 20 DEGs encoding bZIP-family proteins, of which 17 were upregulated, bZIP60 significantly so, suggesting that bZIPs play important roles in sensing ER stress and activating gene expressing in P. tomentosa under heat stress. In addition, MYB27 and MYB4R1 were upregulated after 5 h of heat stress, whereas WRKY50 was upregulated after 12 h in our qRT-PCR validation experiments. Therefore, P. tomentosa may employ MYBs and WRKYs to response heat at different phases after heat stress.

Conclusions
Heat stress induced significant physiological, transcriptomic, and metabolomic changes in the leaves of P. tomentosa. Heat stress increased cellular H 2 O 2 , antioxidant enzyme activities, proline, raffinose, and melibiose, whereas it decreased soluble sugar, pyruvate, fumarate, and myo-inositol. In photosynthesis, heat stress suppressed expressions of most genes related to the light-harvesting complexes and photosynthetic electron transport system. In respiration, heat stress enhanced glycolysis, and retarded the TCA cycle. Heat stress significantly affected carbohydrate metabolic pathways. Our results provide important information of metabolites, genes and transcription factors that involved in the response to heat stress in P. tomentosa.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1999-4907/10/5/383/s1, Figure S1: Volcano Plot for DEGs (Heat vs CTRL), Figure S2: Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. (a) Top 30 enriched GO terms between heat treatments and the control (arrows indicate important GO terms associated with heat stress). * Significant enrichment.
(b) Top 20 enriched KEGG pathways between heat treatments and the control (arrows, significant enrichments associated with heat stress), Table S1: A list of the quality of sequencing data, Table S2: Detailed information of some important DEGs under heat stress for 5 hours, Table S3: Detailed information of DEGs (more than 1.5-log2 fold change) under heat stress, Table S4: DEGs encoding transcription factor families, Table S5: Primers used for qRT-PCR.
Author Contributions: S.R. and B.J. designed the research. K.M., Z.L. and S.R. performed the experiments. J.C., G.C., N.T. and P.T. participated to the result analyses and interpretation of data. S.R., B.J. and L.W. wrote the manuscript.