Transcriptomic analysis reveals unique molecular factors for lipid hydrolysis, secondary cell-walls and oxidative protection associated with thermotolerance in perennial grass

Background Heat stress is the primary abiotic stress limiting growth of cool-season grass species. The objective of this study was to determine molecular factors and metabolic pathways associated with superior heat tolerance in thermal bentgrass (Agrostis scabra) by comparative analysis of transcriptomic profiles with its co-generic heat-sensitive species creeping bentgrass (A. stolonifera). Results Transcriptomic profiling by RNA-seq in both heat-sensitive A. stolonifera (cv. ‘Penncross’) and heat-tolerant A. scabra exposed to heat stress found 1393 (675 up- and 718 down-regulated) and 1508 (777 up- and 731 down-regulated) differentially-expressed genes, respectively. The superior heat tolerance in A. scabra was associated with more up-regulation of genes in oxidative protection, proline biosynthesis, lipid hydrolysis, hemicellulose and lignin biosynthesis, compared to heat-sensitive A. stolonifera. Several transcriptional factors (TFs), such as high mobility group B protein 7 (HMGB7), dehydration-responsive element-binding factor 1a (DREB1a), multiprotein-bridging factor 1c (MBF1c), CCCH-domain containing protein 47 (CCCH47), were also found to be up-regulated in A. scabra under heat stress. Conclusions The unique TFs and genes identified in thermal A. scabra could be potential candidate genes for genetic modification of cultivated grass species for improving heat tolerance, and the associated pathways could contribute to the transcriptional regulation for superior heat tolerance in bentgrass species. Electronic supplementary material The online version of this article (10.1186/s12864-018-4437-z) contains supplementary material, which is available to authorized users.

One approach to unraveling mechanisms of plant tolerance to stresses is to examine plants adapted to extremely stressful environments. A temperate (C3) perennial grass species, thermal bentgrass (A. scabra) endemic to geothermal areas of Yellowstone National Park, exhibits superior heat tolerance to other C3 grass species, as it is able to survive at soil temperature up to 45°C [12,13], while soil temperature over 18°C or air temperature over 24°C is detrimental for most C3 grass species [14]. Physiological, proteomic, and metabolic analysis with thermal bentgrass have found that superior heat tolerance of A. scabra was associated with the adjustment of various metabolic processes, including lowering respiratory consumption of carbohydrates, increases of alternative respiration and carbon use efficiency [15][16][17][18], activation of antioxidant metabolism, induction of stress-protective proteins, such as heat shock proteins [19][20][21] and the accumulation of osmoprotectants, such as soluble sugars and proline [22]. However, the molecular factors underlying the superior heat tolerance of the thermal grass species are not well documented, but such information is useful for improving heat tolerance in cultivated grass species.
The objective of this study was to identify unique transcriptional factors and genes, as well as the associated metabolic pathways accounting for the superior heat tolerance of the wild grass species, thermal A. scabra, by comparative analysis of the transcriptomic changes in response to heat stress between thermal A. scabra and its co-generic heat-sensitive species (A. stolonifera).

Plant materials and growth conditions
Tillers (30 per individual plant) of A. stolonifera ('Penncross') or A. scabra ('NTAS') were collected from stock plants and transferred to plastic containers (57 × 44 × 30 cm, 12 drainage holes) filled with fritted clay medium (Profile Products, Deerfield, IL). Plants were established for 35 d in a greenhouse with average temperature of 23/20°C (day/night), 60% relative humidity (RH), and 750 μmol m −2 s −1 photosynthetically active radiation (PAR) from natural sunlight and supplemental lighting. Plants were irrigated daily, fertilized twice per week with half-strength Hoagland's nutrient solution [23], and trimmed to 2 cm once per week during establishment. Plants were not trimmed during the final week of establishment to allow for sufficient regrowth prior to stress imposition, after which time all plants were transferred to controlled-environment growth chambers (Environmental Growth Chamber, Chagrin Falls, Ohio).

Heat stress treatments and experimental design
Plants were maintained in controlled-environment growth chambers controlled at 22/18°C (day/night), 600 μmol m −2 s −1 PAR, 60% RH, and 14-h photoperiod for one week prior to stress imposition, and then air temperature was raised to 35/30°C to impose heat stress for 21 d. During stress treatment, plants were irrigated daily, and fertilized twice per week with half-strength Hoagland's nutrient solution. The experiment was arranged in a split-plot design with temperature treatment (control or heat) as the main plots and grass species (A. scabra or A. stolonifera) as subplots. Each species was replicated in four containers and each temperature treatment was repeated in four growth chambers. Plants under the same temperature were relocated across growth chambers every 3 d to avoid possible confounding effects of chamber environmental variations.

Physiological measurements
Leaf relative water content (RWC), chlorophyll content (Chl) and electrolyte leakage (EL) were measured at 0 and 21 d of heat stress to assess differential physiological responses of the two plant species under both control and heat stress conditions. Approximately 0.8 g fresh leaf tissue was collected from four individual plants per line per container, and then pooled for RWC, EL, and Chl measurements. For RWC, 0.2 g of leaf blades were first weighed for fresh weight (FW), soaked in water for 12 h and again weighed for turgid weight (TW), dried in an oven at 80°C for 3 d, and finally weighed for dry weight (DW). RWC was calculated using the formula (%) = ([FW -DW] / [TW -DW]) × 100 [24]. For Chl, approximately 0.2 g fresh leaf tissue was submerged in 10 ml dimethyl sulphoxide for 3 d to extract total chlorophyll. The absorbance of the leaf extract was measured at 663 nm and 645 nm with a spectrophotometer (Spectronic Genesys 2; Spectronic Instruments, Rochester, NY) and Chl calculated using the formula described in [25]. For EL, approximately 0.2 g of fresh leaf tissue was rinsed with deionized water, placed in a test tube containing 30 mL deionized water, agitated on a conical shaker for 12 h, and initial conductance (Ci) measured using a conductivity meter (YSI Model 32, Yellow Springs, OH). Tubes containing leaf tissue were then autoclaved at 121°C for 20 min and again agitated for 12 h. The maximal conductance (C max ) of incubation solution was then measured and EL (%) was calculated as ((C i /C max ) × 100) [26]. Four biological replicates (n = 4) of each species were performed for each parameter under either control or heat stress condition, respectively. Statistical differences between treatment means were separated by Student's t-test at the P level of 0.05.
Gene ontology (GO) term classification was performed by CateGOrizer [32], using "GO_slim2" method. The GO enrichment analysis for DEGs was performed using GOEAST [33], by first implementing Customized Result Analysis for up-and down-regulated DEGs in each species, respectively, and then comparing between two species in Multi-GOEAST, using default parameters. KEGG pathway enrichment analysis was performed using DAVID v6.8 [34], by using UniProt IDs for the entire transcriptome background and DEGs in both species.
The transcriptome shotgun assembly of both A. stolonifera and A. scabra were deposited at GenBank Transcriptome Shotgun Assembly (TSA) database, under the accession of GFJH00000000 and GFIW00000000, respectively. The version described in this paper is the first version, GFJH01000000 and GFIW01000000.

Validation of gene expression levels
Gene expression analysis was performed by quantitative reverse transcriptase polymerase chain reaction (qRT-PCR). Total RNA was isolated from ground leaf powder using TRIzol reagent (Life Technologies, Grand Island, NY) and treated with DNase (TURBO DNA-free kit; Life Technologies, Grand Island, NY) to remove contaminating genomic DNA. Total RNA (2 μg) was reverse-transcribed using a high-capacity cDNA reverse transcription kit (Life Technologies, Grand Island, NY). The synthesized cDNA was amplified in a StepOnePlus Real-Time PCR system (Life Technologies, Grand Island, NY) using the following parameters: pre-heat cycle of 95°C for 3 min, 40 cycles of 95°C denaturation for 30 s per cycle, and 60°C annealing/extension for 30 s per cycle. Power SYBR Green PCR Master Mix (Life Technologies, Grand Island, NY) was the intercalating dye used to detect gene expression level. Gene name, accession number, forward and reverse primer sequences are provided in Table 1. A melting curve analysis was performed for each primer set to confirm its specificity. Actin was used as the reference gene, since its expression was consistent throughout treatments. A ΔΔCt method was used to calculate the relative expression level between genes of interest and reference gene, respectively [35]. Four biological replicates (n = 4) from each species were performed for each gene under either control or heat stress condition, respectively. Statistical differences between treatment means were separated by Student's t-test at the P level of 0.05.

Physiological responses to heat stress
Under control conditions, leaf relative water content (RWC) did not differ significantly between A. stolonifera and A. scabra. Heat treatment caused significantly decline in RWC at 21 d in both A. stolonifera and A. scabra, by 13.1% and 5.8%, respectively. However, RWC in A. scabra was significantly higher than that in A. stolonifera (Fig. 1a). No significant differences in leaf chlorophyll content (Chl) were found between A. stolonifera and A. scabra under control conditions. At 21 d of heat treatment, Chl content decreased significantly in both A. stolonifera and A. scabra, by 17.8% and 9.6%, respectively; leaf Chl in A. scabra was significantly higher than that in A. stolonifera (Fig. 1b). For electrolyte leakage (EL), there was no significant difference found between A. stolonifera and A. scabra under control conditions. Heat stress at 21 d resulted in significantly increases in EL in both A. stolonifera and A. scabra, by 63.7% and 47.6%, respectively. Leaf EL in A. scabra was significantly lower than that in A. stolonifera (Fig. 1c).
Next-generation sequencing of A. stolonifera and A. scabra The RNA sequencing yielded more than 19 million reads per library of A. stolonifera and A. scabra plants exposed to non-stress control and heat stress conditions, providing over 5× coverage of the estimated genome of A. stolonifera ( Table 2). The de novo transcript assembly by Trinity algorithm had good alignment rate, indicating that the assembled transcripts were largely representing transcriptome in these two species (Table 2). In addition, transcript qualities were also confirmed by long N50 numbers, contig lengths and similar GC contents (Table 3). It is therefore indicated that the Illumina RNA-seq was successfully performed to obtain transcriptional profiles for A. stolonifera and A. scabra under heat stress.
After transcript clustering and annotation, a total of 75,253 and 81,597 UniGenes were obtained by BlastX against NCBI protein NR database (Table 4). Further annotation with GO, KEGG, COG and Pfam also had similar results among them. The components of annotation were mainly from Arabidopsis and rice (Table 5). GO term classification showed that the functional distributions of UniGenes were similar between A. stolonifera and A. scabra (Fig. 2).

GO term enrichment analysis
Using the threshold of |log2FC| > 1, and FDR < 0.01, we identified 675 and 777 up-regulated DEGs, and 718 and 731 down-regulated DEGs in A. stolonifera    and A. scabra, respectively, by heat stress (Fig. 3). In order to find out specific molecular factors for the superior heat tolerance in A. scabra, up-and downregulated DEGs due to heat tress were analyzed by GO term enrichment analysis in the two species separately (Figs. 4, 5, 6, 7, 8 and 9; For heat map, see Additional files 1 and 2). In the up-regulated DEGs, several functional categories were enriched only in A. scabra, including hemicellulose metabolic process, cell wall biogenesis, L-proline biosynthetic process, lipid catabolic process, lipid transport, lignan biosynthetic process for Biological Process (BP) terms (Fig. 4); In Molecular Function (MF) terms, monooxygenase activity, oxidoreductase activity, several glucosidase activity, and several monosaccharidase activity, such as arabinosidase activity, mannosidase activity, galactosidase activity, fucosidase activity were also uniquely enriched in A. scabra (Fig. 5). The uniquely enriched DEGs of A. scabra in Cellular Component (CC) terms were mainly at anchored component of membrane and apoplast region (Fig. 6). Some down-regulated DEGs were found to be enriched only in A. stolonifera, including DNA-templated transcription, glucose metabolic process, several amino acid metabolic process, such as L-serine, cysteine, and glycine, pentose-phosphate shunt, hydrogen peroxide catabolic process, chloroplast organization, regulation of photosynthesis, positive regulation of translation, and response to oxidative stress in BP terms (Fig. 7). Several cofactor binding functions, such as poly(U) binding, NAD binding, NADP binding, FMN binding, beta-glucosidase activity, cis-trans isomerase activity, several transaminase activity, sulfate adenyltransferase (ATP) activity, adenylate kinase activity, transketolase activity, glyceraldehyde-3-phosphate dehydrogenase (GAPDH) activity, glycolate oxidase activity, glucose-6-phosphate dehydrogenase activity, monooxygenase activity and peroxidase activity were also uniquely enriched in A. stolonifera in MF terms (Fig. 8). The CC terms further showed that down-regulated transcripts uniquely enriched in A. stolonifera were located in oxidoreductase complex, apoplast, NAD(P)H dehydrogenase complex, peroxisome, and chloroplast membrane (Fig. 9). The biological process and molecular functions of GO terms in up-regulated DEGs showing specific enrichment to A. scabra, and the GO terms in down-regulated DEGs showing specific enrichment to A. stolonifera were identified, and the individual transcripts in each category were also analyzed (Tables 6 and 7). In the upregulated DEGs, those related to cell wall biogenesis, lipid metabolism, proline biosynthesis, lignan biosynthesis, oxidoreductase activity and glucosidase activity, were uniquely enriched in A. scabra ( Table 6). The down-regulated DEGs found only in A. stolonifera included dhurrin biosynthetic process, amino acid metabolism, glucose metabolic process, pentose phosphate shunt, peroxidase activity, beta-glucosidase activity, cistrans isomerase activity, aminotransferase activity, sulfate adenylyltranferase activity, transketolase activity, Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) activity, glycolate oxidase activity, chloroplast organization, regulation of transcription and translation, energy metabolism and monooxygenase activity (Table 7).

KEGG pathway enrichment analysis
KEGG pathway enrichment analysis compared DEGs between A. stolonifera and A. scabra, and identified pathways in the degree of enrichment upon heat stress (Tables 8 and 9). In the up-regulated DEGs by heat stress, the top six enriched pathways in A. scabra were cutin, suberine and wax biosynthesis, biosynthesis of secondary metabolites, metabolic pathways, fatty acid elongation, phenylpropanoid biosynthesis, ABC transporters, and those in A. stolonifera were biosynthesis of secondary metabolites, arginine and proline metabolism, alpha-linolenic acid metabolism, galactose metabolism, beta-alanine metabolism and plant-pathogen interaction ( Table 8). In the down-regulated DEGs, the top six enriched pathways were the same in both A. stolonifera and A. scabra, including metabolic pathways, biosynthesis of secondary metabolites, glyoxylate and dicarboxylate metabolism, carbon metabolism, glycine, serine and threonine metabolism and biosynthesis of antibiotics (Table 9). Biological Process (BP) of GO term enrichment for up-regulated DEGs in A. stolonifera and A. scabra. Green color indicates GO terms that were specifically enriched in A. stolonifera. Red color indicates GO terms that were specifically enriched in A. scabra. Yellow color indicates GO terms that were commonly enriched in both A. stolonifera and A. scabra. The density of color was proportional to statistical significance, which was shown as p1 for P-value of A. stolonifera and p2 for P-value of A. scabra Green color indicates GO terms that were specifically enriched in A. stolonifera. Red color indicates GO terms that were specifically enriched in A. scabra. Yellow color indicates GO terms that were commonly enriched in both A. stolonifera and A. scabra. The density of color was proportional to statistical significance, which was shown as p1 for P-value of A. stolonifera and p2 for P-value of A. scabra

Transcription factors related to heat tolerance
Transcription factors (TFs) responsive to heat stress showed high similarity between A. stolonifera and A. scabra, including up-regulation of ABA-inducible, basic Helix-Loop-Helix (bHLH), ethylene-responsive factor (ERF), protein FD, G-box-binding factor, heat stress factor (HSF), homeobox-leucine zipper, MYB, NAC, nuclear transcription factor Y, WRKY, and down-regulation of Green color indicates GO terms that were specifically enriched in A. stolonifera. Red color indicates GO terms that were specifically enriched in A. scabra. Yellow color indicates GO terms that were commonly enriched in both A. stolonifera and A. scabra. The density of color was proportional to statistical significance, which was shown as p1 for P-value of A. stolonifera and p2 for P-value of A. scabra APG, PHL1-like, RNA-polymerase sigma factor, zincfinger protein (Table 10). However, some TFs were uniquely regulated by heat stress in A. scabra, such as the up-regulation of high mobility group B protein 7 (HMGB7), dehydration-responsive element-binding factor 1A (DREB1A), multiprotein-bridging factor 1c, CCCHdomain containing protein 47, and down-regulation of GLK1, GATA transcription factor 21 and 26, protein REVEILLE, ASR3, HY5 (Table 11).

Validation of RNA-seq with qRT-PCR
The differential expressions of several DEGs in the RNA-seq data were validated using qRT-PCR. Heat stress significantly increased gene expression levels of xyloglucan endo-transferase 25 (XET25), GDSL esterase, Dirigent protein 5, pyrroline-5-carboxylate reductase (P5CR), Cytochrome P450 77A3, HMGB7, and DREB1A in both A. stolonifera and A. scabra. However, expression levels for these genes in A. scabra under heat stress were significantly higher than those in A. stolonifera ( Fig. 10a-g). Heat stress significantly decreased gene expression levels of glycine cleavage system H protein, GAPDH A, peroxidase 4, and beta-glucosidase 3 in both A. stolonifera and A. scabra. However, the expression levels for these genes in A. scabra under heat stress were significantly higher than those in A. stolonifera (Fig. 10h-k). Heat stress also significantly increased expression level of transcription factor DIVARICATA in both A. stolonifera and A. scabra, but the expression level in A. stolonifera under heat stress was significantly higher than that in A. scabra (Fig. 10l). Furthermore, the fold changes of these genes obtained by qRT-PCR analysis were compared with RNA-seq data. The Person's correlation coefficient between data from RNA-seq and qRT-PCR was 0.95 for A. scabra, and 0.93 for A. stolonifera, indicating that the transcriptional regulations under heat stress in these two species were valid regardless of detecting methods (Table 12).

Discussion
The comparative analysis of transcriptome profiles between A. stolonifera and A. scabra exposed to heat stress found that metabolic processes involved in heat responses were similar in the two species, but some TFs and genes uniquely enriched in A. scabra, which could account for its superior heat tolerance. The following sections focus on the discussion of uniquely upregulated TFs and genes in A. scabra and uniquely down-regulated TFs and genes in A. stolonifera regarding their functions and roles in heat tolerance.
Previous studies with A. scabra found that accumulation of carbohydrates and amino acids play major roles in heat tolerance for this heat-tolerant grass Fig. 7 Biological Process (BP) of GO term enrichment for down-regulated DEGs in A. stolonifera and A. scabra. Green color indicates GO terms that were specifically enriched in A. stolonifera. Red color indicates GO terms that were specifically enriched in A. scabra. Yellow color indicates GO terms that were commonly enriched in both A. stolonifera and A. scabra. The density of color was proportional to statistical significance, which was shown as p1 for P-value of A. stolonifera and p2 for P-value of A. scabra Green color indicates GO terms that were specifically enriched in A. stolonifera. Red color indicates GO terms that were specifically enriched in A. scabra. Yellow color indicates GO terms that were commonly enriched in both A. stolonifera and A. scabra. The density of color was proportional to statistical significance, which was shown as p1 for P-value of A. stolonifera and p2 for P-value of A. scabra species [20,22,36,37]. The GO term and KEGG analysis identified some specific pathways of genes involved in carbohydrate, amino acid, and energy metabolism, including unique up-regulation of glucosidases and monosaccharidases activity in A. scabra (Fig. 5, Table 6), and down-regulation of aminotransferase activity, glucose metabolic process, pentose phosphate shunt, transketolase, cis-trans isomerase, GAPDH, glycolate oxidase activity, cofactor binding (NAD, NADP, FMN), glucose-6phosphate dehydrogenase and chloroplast organization in A. stolonifera (Fig. 8, Table 7). In addition, serine hydroxymethyltransferase 1 (SHMT1) was significantly downregulated only in A. stolonifera (Table 7). SHMT catalyzes the interconversion between serine and glycine [38]. Previous study of root proteomic profiles between A. stolonifera and A. scabra under heat stress showed that one SHMT protein spot was decreased only in A. stolonifera, which agreed our transcriptional observation [19,20]. In contrast to transcript responses of heatsensitive A. stolonifera, lack of down-regulation of some genes mentioned above in A. scabra suggest that the maintenance of transcriptional levels of genes in carbohydrate, and amino acid, and energy metabolism may be associated with the corresponding metabolite accumulation under heat stress, contributing to the superior heat tolerance.
Under heat stress, A. scabra showed up-regulation of several functional categories that were related in antioxidative responses and antioxidant protection, while many of the functional categories related to oxidative protection, such as peroxidase activity, peroxisome, were down-regulated in A. stolonifera (Figs. 8 and 9). Most of the up-regulated antioxidant-related genes were Cytochrome P450s ( Table 6). The cytochrome P450 is a superfamily catalyzing various oxidative reactions, including biosynthesis of lipophilic compounds (fatty acids, sterols, cutin, suberine and wax, phenylpropanoids, brassinolides and gibberellins [39]. The microarray analysis of cytochrome P450 family in Arabidopsis showed that they are highly responsive to both abiotic and biotic stresses [40]. Other up-regulated genes involved in antioxidant defense included oxidoreductase and monooxygenase activity found in A. scabra under heat stress, which were also mainly involved in plant antioxidative response. ROS content in A. scabra root tissue under heat stress was significantly lower than that in A. stolonifera [41], suggesting that the ROS scavenging capacity was better maintained in A. scabra under heat stress. It is also worthy to point out the 4hydroxyphenylacetaldehyde oxime monooxygenase that were uniquely down-regulated in A. stolonifera was also a Cytochrom P450 gene (CYP71E1), which is involved in Fig. 9 Cellular Component (CC) of GO term enrichment for down-regulated DEGs in A. stolonifera and A. scabra. Green color indicates GO terms that were specifically enriched in A. stolonifera. Red color indicates GO terms that were specifically enriched in A. scabra. Yellow color indicates GO terms that were commonly enriched in both A. stolonifera and A. scabra. The density of color was proportional to statistical significance, which was shown as p1 for P-value of A. stolonifera and p2 for P-value of A. scabra  the oxime-metabolizing step in biosynthesis of dhurrin [42]. Little information was known about dhurrin and its relation to heat response in plants, which deserves further investigation. Transcripts in proline biosynthesis, mainly pyrroline-5-carboxylate reductase (P5CR), were also up-regulated in A. scabra under heat stress (Fig. 4, Table 6). P5CR is the final step in proline biosynthesis pathway, which reduces proline-5-carboxylate to proline [43]. It is generally accepted that proline acts as a cellular osmolyte, and thus its accelerated biosynthesis indicates enhanced plant osmotic stress resistance [44]. In addition, proline is also involved in maintenance of redox balance and ROS scavenging [45][46][47]. Higher levels of proline were identified in heat-stressed cotton (Gossypium hirsutum L.) [48], and its positive role in heat tolerance was confirmed in various plant species [1,49,50]. Our previous study found that proline content was significantly higher in A. scabra root tissues under heat stress than that in A. stolonifera [22]. However, there is little information regarding to P5CR expression regulation under heat stress. This study found the up-regulation of pyrroline-5-carboxylate reductase in A. scabra, which could play positive roles in the maintenance of proline synthesis in the heat-tolerant species under heat stress.
Most of the transcripts up-regulated in lipid catabolic process were GDSL esterases and Phospholipase A1 ( Table 6). The GDSL-motif enzyme is a newly discovered lipase family that shares the highly conserved motif Gly-X-Ser-X-Gly (X means any amino acid) in the sequence [51,52]. The number of GDSL esterase/lipase family members ranged from 57 to 130 in several plant organisms [53,54]. The GDSL esterases/lipases might play an important role in plant development and morphogenesis [52]. Some of the GDSL esterases were reported to confer plant abiotic stress tolerance, such as drought and salt stress [55,56]. Phospholipase A1 is one of the multigene family of phospholipases, hydrolyzing the sn-1 acylester bond of phospholipids to free fatty acids and 2-acyl-1-lysophospholipids [57]. Compared to mammalian phospholipase A1s, only a few genes were discovered in plants. A phospholipase A1 homolog in Arabidopsis,      The transcriptional regulations under heat stress, log2 fold change (log2 FC), in these GO terms are also listed AtDAD1, was placed in the initial step of jasmonic acid biosynthesis, making it important for plant responses to abiotic stress, tendril coiling, fruit ripening and developmental maturation of stamens and pollens [58]. Another phospholipase A1 in hot pepper (Capsicum annuum) showed high sequence similarity to Arabidopsis [59]. Another phospholipase A1 homolog in Arabidopsis, AtLCAT3, was determined to have in vitro enzymatic activity, although its molecular function has yet to be assigned [60]. Therefore, the GDSL esterases and phospholipase A1s found in the up-regulated transcripts in A. scabra were also considered to be involved in lipid catalysis, possibly through jasmonic acid signal transduction pathway. This is first report of GDSL esterases and phospholipase A1s related to heat tolerance. Further studies regarding to their functions and regulation of heat tolerance in plants are needed.
The transcripts involved in cell wall structure and properties were up-regulated in A. scabra under heat stress, including xyloglucan endo-transglycosylases (XETs), and cellulose synthase (Table 6). XETs make nonhydrolytic cleavage and ligation of xyloglucan chains, which is involved in cell wall loosening [61]. Cellulose synthase family is also well-defined, and involved in the formation of plant primary and secondary cell wall [62] Plant cell wall structure undergoes reassembly that involves biosynthesis of major cell wall components during plant responses to abiotic stress [63][64][65]. Xu et al. [66] reported that transcript levels of XETs in tall fescue root tissues were decreased under water stress, and exogenous application of ascorbic acid could mitigate the reduction. Little information was known regarding to genes for cell wall biosynthesis and properties related to heat tolerance.
Several transcripts involved in lignan biosynthetic process were also up-regulated in A. scabra, including dirigent protein 5, aureusidin synthase 1, and      (Table 6). Dirigent proteins play an important role in monolignol coupling to both lignin and lignan formations [67]. Dirigent protein family was reported to participate in defense responses, secondary metabolism, temperature, and salinity stress [68][69][70]. Aureusidin synthase is a binuclear copper enzyme, and a homolog of plant phenol oxidase [71]. It is proposed to be a chalcone-specific plant phenol oxidase for aurone biosynthesis [72]. Larreatricin hydroxylase is an enantio-specific polyphenol oxidase [73], but their physiological roles in plant adaptation to abiotic stress are unknown. Results in our study indicated that the up-regulation of genes involved in secondary cell-wall materials could contribute to the maintenance of cell wall structure and functional properties for A. scabra to maintain growth under heat stress. Several transcription factors were uniquely up-regulated under heat stress, including high mobility group B protein (HMGB) 7, dehydration-responsive element-binding factor (DREB) 1a, Multiprotein-bridging factor (MBF) 1c, and CCCH-type zinc finger protein 47 (Table 11). The high mobility group B protein (HMGB) belongs to chromatin-associated proteins, and acts primarily as architectural facilitator in nucleoprotein complex assembly and transcriptional regulation and recombination [74,75]. Little is known about its function in plant stress responses, except that Arabidopsis HMGBs showed induced expression levels under cold stress [76]. Dehydration-responsive element-binding factor (DREB) is one of the sub-groups in AP2/EREBP family, and activates target genes that have dehydrationresponsive elements (DREs) [77,78]. DREB1 and DREB2 were reported to confer plant drought, salinity and low-temperature tolerance [79][80][81], while only Arabidopsis DREB2A has dual functions in water-stress and heat-stress response [82]. Multiprotein-bridging factor 1 (MBF1) is a conserved transcriptional coactivator that bridges a basic region/leucine zipper (bZIP) type coactivator and a TATA-box binding protein [83,84]. The Arabidopsis MBF1c was induced in response to heat stress [84,85]. The transgenic Arabidopsis constitutively expressing MBF1c enhanced plant heat tolerance by perturbing ethylene response signal transduction pathway [86]. Plant CCCH-type zinc finger proteins were shown to be involved in embryo formation, floral reproductive organ formation, delay of leaf senescence, and calmodulin-mediated RNA processing [87][88][89][90]. Two CCCH-type zinc finger proteins, AtSZF1 and AtSZF2, were induced upon salt stress, and negatively regulate salt-responsive genes in Arabidopsis [91]. Our results suggest that the up-regulation of HMGB7, DREB1a, MBF1c, and CCCH-type zinc finger protein 47 could active the related down-stream genes regulating heat tolerance. Further research is needed to identify the down-stream genes in order to unravel the molecular roles of those transcriptional factors in heat tolerance.

Conclusions
In summary, our comparative analysis of transcriptomic changes in response to heat stress for heat-tolerant thermal A. scabra and heat-sensitive A. stolonifera showed divergent transcriptional regulations of heat tolerance in perennial grass species, which complemented to previous findings of physiological traits and proteins conferring the superior heat tolerance of the thermal species adapted to extremely high soil temperature. The