Metabolite and transcript profiling of Guinea grass (Panicum maximum Jacq) response to elevated [CO2] and temperature

Introduction By mid-century, global atmospheric carbon dioxide concentration ([CO2]) is predicted to reach 600 μmol mol−1 with global temperatures rising by 2 °C. Rising [CO2] and temperature will alter the growth and productivity of major food and forage crops across the globe. Although the impact is expected to be greatest in tropical regions, the impact of climate-change has been poorly studied in those regions. Objectives This experiment aimed to understand the effects of elevated [CO2] (600 μmol mol−1) and warming (+ 2 °C), singly and in combination, on Panicum maximum Jacq. (Guinea grass) metabolite and transcript profiles. Methods We created a de novo assembly of the Panicum maximum transcriptome. Leaf samples were taken at two time points in the Guinea grass growing season to analyze transcriptional and metabolite profiles in plants grown at ambient and elevated [CO2] and temperature, and statistical analyses were used to integrate the data. Results Elevated temperature altered the content of amino acids and secondary metabolites. The transcriptome of Guinea grass shows a clear time point separations, with the changes in the elevated temperature and [CO2] combination plots. Conclusion Field transcriptomics and metabolomics revealed that elevated temperature and [CO2] result in alterations in transcript and metabolite profiles associated with environmental response, secondary metabolism and stomatal function. These metabolic responses are consistent with greater growth and leaf area production under elevated temperature and [CO2]. These results show that tropical C4 grasslands may have unpredicted responses to global climate change, and that warming during a cool growing season enhances growth and alleviates stress. Electronic supplementary material The online version of this article (10.1007/s11306-019-1511-8) contains supplementary material, which is available to authorized users.


Introduction
Atmospheric carbon dioxide concentrations ([CO 2 ]) have risen by 40% since pre-industrial times, with approximately half of that increase within the past 40 years. The International Panel on Climate Change (IPCC) has predicted that [CO 2 ] will rise to 600 ppm and mean global temperature will increase 1.1 to 2.6 °C by the end of the century, with a moderate climate change projection scenario (IPCC 2014). Although the combined effects of rising atmospheric [CO 2 ] and rising temperatures have been studied in field experiments in temperate forests Svensson et al. 2018), croplands (Ruiz-Vera et al. 2015;Cai et al. 2016), and grasslands (Ryan et al. 2015), there have been fewer field studies of the combined effects of rising [CO 2 ] and temperature in tropical and sub-tropical ecosystems. Differences in climate and nutrient limitation in tropical regions may fundamentally alter how tropical ecosystems respond to global climate change compared to well-studied temperate ecosystems (Leakey et al. 2012). Therefore, the paucity of data from tropical regions presents a significant challenge for an accurate understanding of global responses to climate change. Adding to this challenge is the realization that responses of vegetation to combined elevated [CO 2 ] and temperature treatments are different from single-factor experiments, and that additive effects are rare (Dieleman et al. 2012). Thus, there is a need for multi-factor global change experiments, especially in tropical regions.
Pastures cover about 30% of the terrestrial surface in tropical regions and play an important role in carbon sequestration (Boval and Dixon 2012;Ramankutty et al. 2008). In tropical regions, C 4 grasses dominate those pastures (Still et al. 2003). Panicum maximum Jacq. (Guinea grass) is a perennial C 4 grass species native to tropical Africa, which has been adapted in Brazil as a forage and range grass (Pedreira et al. 2015). Approximately 90% of Brazilian cattle are raised on pasture land (Cardoso et al. 2016), and perennial C 4 grasses such as Guinea grass are principle species for the Brazilian feedstock industry. Despite the economic importance of Guinea grass to this industry, very little is known about how it and other C 4 forage grasses will respond to rising atmospheric [CO 2 ] and temperature. To address this knowledge gap, a Temperature Free Air CO 2 Enrichment (T-FACE) experiment was established in Ribeirão Preto, Brazil to expose Guinea grass to elevated [CO 2 ] (600 ppm) and increased canopy temperatures of 2 °C (Britto de Assis Prado et al. 2016).
'Omic' technologies including transcriptomics and metabolomics provide a functional analysis to connect physiological and molecular responses to genetic and phenotypic information (Langridge and Fleury 2011). Such technologies can be used to better understand the mechanisms underpinning plant responses to global climate change and can reveal targets for improving plant performance in the future. In recent years there has been an increasing number of studies investigating transcriptomic and metabolomic responses of C 4 grasses to a variety of abiotic stresses imposed in greenhouses and/or chambers (e.g., Sicher et al. 2012;Meyer et al. 2014;Sui et al. 2015;Toledo-Silva et al. 2013). However, there is an increasing need to investigate and understand plant transcriptional and metabolic responses to global change or abiotic stress in the field where responses can be dramatically different from controlled environments (Lovell et al. 2016). The combined use of FACE and 'omic' technologies provides a powerful approach for gaining new insight into plant molecular responses to global climate change in the field. Combining metabolomic and transcriptomic studies allowing for a holistic view of the plant's response to stress.
Our study combines the power of FACE experiments for growing plants under elevated [CO 2 ] and temperature in fully open-air conditions with transcriptomic and metabolomics profiling. The aims are to (1) understand the biochemical responses of Guinea grass to rising [CO 2 ] and temperature when applied individually and in combination; (2) understand how gene expression is altered in each treatment; and (3) integrate the metabolomic and transcriptional responses.

Experimental design and tissue collection
The field experiment was conducted in a 2500 m 2 field at the University of São Paulo in Ribeirão Preto, São Paulo Brazil. Guinea grass was planted in a completely replicated design with 4 ambient temperature and [CO 2 ] plots (C), 4 elevated [CO 2 ] (600 μmol mol −1 ) and ambient temperature plots (eC), 4 ambient [CO 2 ] and elevated temperature [+ 2 °C] plots (eT), and 4 combined treatment (elevated [CO 2 ] and elevated temperature) plots (eCeT). Each plot was 2-m in diameter; [CO 2 ] was controlled using Free Air CO 2 Enrichment, and canopy warming was provided by infrared ceramic heaters as described by Britto de Assis Prado et al. (2016). Treatment was started on April 22, 2014.
Leaf tissue was collected for transcriptomic and metabolomics analysis on two dates: May 22, 2014 (time point A-30 days of treatment), and June 14, 2014 (time point B-50 days of treatment). Six leaves per experimental plot were collected at each time point and then pooled for RNA extraction and metabolomics analysis. At the time of sampling, tissue was immediately quenched in liquid nitrogen. RNA was extracted in Brazil using the protocol described by Bilgin and colleagues (Bilgin et al. 2009), and subsequently shipped in ethanol to the U.S. for transcriptomic analysis. Lyophilized leaf tissue from the same leaves was sent for metabolite analysis. Leaves were sampled for biomass, area, and specific leaf area on August 22 and 29th, 2014 from 60 tillers per treatment as described by Britto de Assis Prado et al. (2016).

Transcriptional analysis
RNA samples that were shipped from Brazil were resuspended by washing with sodium acetate (pH 5.5) and glycerol. RNA quantity and quality were re-assessed using a spectrophotometer (Nanodrop 1000, Thermo Fisher) and microfluidic visualization tool (Bioanalyzer, Agilent Technologies). 32 samples (16 per time point including 4 from each [CO 2 ] and temperature treatment) were submitted for library construction and sequencing at the Roy J. Carver Biotechnology Center at the University of Illinois, Urbana-Champaign. RNAseq libraries were prepared with Illumina's TruSeq Standard RNA sample prep kit according to the manufacturer's instruction (Illumina, San Diego, CA, USA). No public reference genome for Guinea grass was available at the time of the sequencing. Therefore, a combination of MiSeq and HiSeq was used to generate sequences of different lengths for de novo genome assembly.
For MiSeq analysis, 50 µg of total RNA per sample was pooled for library preparation without RNA fragmentation. The MiSeq library was quantified by qPCR and sequenced on one MiSeq flowcell for 301 cycles using paired-end sequencing. HiSeq paired-end sequencing was done with four quantified libraries per treatment which were pooled in equimolar concentration and sequenced on two lanes for 161 cycles. The final read lengths for MiSeq and HiSeq were 300 nt and 160 nt in length. A total of 635,649,277 reads were assembled from the MiSeq/HiSeq pools. Quality control for reads generated from sequencing was performed using FastQC (http://www.bioin forma tics.babra ham.ac.uk/ proje cts/fastq c/). Quality reads were used to perform de novo transcriptome assembly using Trinity. The initial assembly consisted of 187,216 genes. A filter was applied to keep only those genes that had at least 10 reads (across the 4 replicates) for an individual treatment. The resulting transcriptome contained 45,073 genes and reads. Functional annotation of the genes was done by using BLAST against Arabidopsis thaliana, Zea mays, and Setaria italica.
Differential gene expression analysis was performed using the R package, LIMMA (Ritchie et al. 2015). Within each time point, individual experimental treatments were compared to the control (ambient [CO 2 ], ambient temperature). Following post 'voom' normalization, extreme outliers were removed, and 39,208 genes were assessed for differential expression analysis. Final normalized reads were log2 transformed for differential gene expression analysis. DESeq 2 (Love et al. 2014) was also run for comparison with the LIMMA approach and results were found to be similar but less conservative (Supplemental Fig. 1). A classical multidimensional scaling plot (MDS) was created using R base package with normalized and transformed reads. MDS plots were used to probe the relationships among samples in multidimensional space using a dissimilarity measure of each pairwise comparison.

Metabolic profiling by GC-MS
Untargeted metabolite analysis was performed with preweighed lyophilized tissue (12.67-19.93 mg, average 16.93 ± 1.92 mg) according to protocol described in Ulanov and Widholm (Ulanov and Widholm 2010). For quality control (QC) 10 µl of leaf extract was taken from each sample and pooled, then run and analyzed after every 9 biological samples. Samples and QC were dried under vacuum and derivatized with 75 µl methoxyamine hydrochloride (Sigma-Aldrich, MO, USA) (40 mg ml −1 in pyridine) for 90 min at 50 °C, then with 125 μl MSTFA + 1%TMCS (Thermo, MA, USA) at 50 °C for 120 min followed by an additional 2-h incubation at room temperature. An internal standard (30 µL hentriacontanoic acid) was added to each sample prior to derivatization. Samples were analyzed on a gas chromatography/mass spectroscopy (GC/MS) system (Agilent Inc, Palo Alto, CA, USA) consisting of an Agilent 7890 gas chromatograph, an Agilent 5975 mass selective detector, and a HP 7683B autosampler. Gas chromatography was performed on a ZB-5MS capillary column (Phenomenex, Torrance, CA, USA). The inlet and MS interface temperatures were 250 °C, and the ion source temperature was adjusted to 230 °C. An aliquot of 1 µl was injected with the split ratio of 10:1. The helium carrier gas constant flow rate was 2.4 ml min −1 . The temperature program was: 5 min isothermal heating at 70 °C, followed by an oven temperature increase of 50 °C min −1 to 310 °C, and a final 10 min at 310 °C. The mass spectrometer was operated in a positive electron impact mode at 69.9 eV ionization energy in m/z 50-800 scan range. The spectra of chromatogram peaks were compared with electron impact mass spectrum libraries NIST08 (NIST, MD, USA), WILEY08 (Palisade Corporation, NY, USA), and a custom library. To allow comparison among samples, data were normalized to the internal standard (hentriacontanoic acid, SIGMA, USA) and standardized to the dry weight of the individual sample. The chromatograms and mass spectra were evaluated using the AMDIS 2.71 (NIST, Gaithersburg, MD, USA) program, using a custom-built database. All known artificial peaks were identified and removed with prior data mining. The instrument variability was within the standard acceptance limit of 5%.
194 metabolites were identified with GC/MS. Technical replicates were averaged across samples. Normalized data was log10 transformed and univariate analysis was performed using SAS (SAS/STAT v9.4, SAS Institute, Inc.) to identify outliers. The MDS plot was created using R base package with normalized and transformed reads. The transformed and normalized metabolomics dataset was analyzed using a general linear model (ANOVA) with SAS PROC GLM where time point was determined non-significant (p > 0.01) and dropped from the model. Means separation tests were done using a two-tailed Dunnett, comparing each treatment to the control. Metabolites were considered significantly different from control if p < 0.01 (Table 1).

Weighted gene correlations network analysis (WGCNA)
WGCNA was run separately for each time point (Langfelder and Horvath 2008

Integration of metabolite and transcriptional data
Linear correlation analysis was performed in R, with the significantly differentially expressed transcripts and metabolites from the second time point (B), along with the physiological data collected and previously reported by Britto de Assis Prado et al. (2016). Linear correlation was selected to highlight possible relationships between the metabolite and transcript abundance. Additionally, the approach taken in time point A, using WGCNA, was not appropriate for the second time point because of high variation within replicates of a given treatment and because far fewer transcripts were significantly different among treatments (see results). Correlations with an adjusted p value (False Discovery rate) of 0.05 or less, and a correlation coefficient of 0.6 or greater were considered significant in this analysis.

Transcriptome analysis
Sequencing produced a total of 635,649,777 paired-end reads, which were initially assembled into 187,216 contigs. Using a filter that kept only transcripts with at least 10 reads across the 4 replicates resulted in a final transcriptome containing 45,073 contigs with an average read length of 400 bp. Annotation of those contigs identified 27,502 unigenes with 11,538 BLASTX hits in the NCBI database. The transcriptome was classified into GO slim terms describing cellular components, molecular functions and biological processes (Supplemental Fig. 2). The highest abundance terms for cellular components, molecular functions, and biological processes were "cell", "catalytic activity", and "metabolic process" respectively, which was similar to another Panicum maximum transcriptome previously published (Toledo-Silva et al. 2013).
Multidimensional scaling (MDS) analysis of leaf transcripts showed a clear time point separation and clustering of treatments within each time point (Fig. 1a). The treatment separation was less apparent in time point B and the variability within the control and eC plots was also greater (Fig. 1a). Differential transcript expression from time point A identified 666 transcripts that were altered in the eCeT treatment when compared to control (p < 0.1, Supplemental File 2). No transcripts were differentially expressed in the eC treatment, when compared to control, and only 3 transcripts were significantly altered in the eT treatment. GO slim terms were assigned to the 666 transcripts altered in the eCeT treatment. Overrepresentation tests were run to identify GO categories with increased occurrence within the differentially expressed genes, with a false discovery rate of 0.05. The most overrepresented GO categories were cellular component "ribosome", molecular function "structural constituent of ribosome", and biological process "protein folding" (Supplemental Fig. 3). Time point B showed the same trend as time point A, but with fewer transcripts significantly changing; only 71 total in eCeT, 2 in eT, and 0 in eC (p < 0.1, Supplemental File 2). The relaxed p value for DEG allowed a larger set of transcripts to test for correlation with metabolite content and/or leaf physiological traits. No GO categories were overrepresented within the significant transcripts in time point B.

Metabolite analysis
Untargeted metabolomics was used to investigate biochemical changes in Guinea grass exposed to eC, eT and eCeT. Although 194 total metabolites were identified, only 125 metabolites were consistently present in the different treatments and those were used for analysis. MDS analysis revealed a separation of the eT treatment (red symbols) from the C plots (black symbols), but no clear separation between sampling dates (Fig. 1b). Additionally, the time points were not significantly different in the ANOVA (p > 0.01). A general linear model with a Dunnett's test was run to test for significant differences in abundance among treatments, and 34 metabolites were significantly altered in at least one treatment compared to the C plots (p ≤ 0.01, Table 1). Only 1-benzylglucopyranoside and isoleucine were significantly different in each of the treatments compared to the C plots.
Overall the metabolome showed the greatest change in the eT treatment compared to C with 30 metabolites significantly changing ( Table 1). All of the metabolites increased in abundance in eT relative to ambient, with the exception of 2-methylmalic acid and oxalic acid (Table 1). Melibiose showed the largest increase in eT amongst sugars, with a 1.94-fold change; arabinose, fructose, galactose, glucose maltose, ribose, and xylose were also significantly increased in eT. Notably sucrose was not altered under any treatment. Many amino acids showed significant increases under eT compared to ambient with Val Thr, and Phe increasing the most (1.64, 1.38, 1.27-fold; Table 1). Organic acids, including neochlorogenic acid and sinapic acid, which are involved in secondary metabolism and specifically monolignol metabolism showed significant increases under the eT treatment compared to C ( Table 1).
The eC treatment did not affect the metabolome as much as eT, with only 2 metabolites showing a log2fold ≥ |1| and p ≤ 0.01. 2-Methylmalic acid (citramalate) had a large decrease (− 3.02-fold) in the eC plots, and this compound decreased across all treatments. Additionally, α-ketoglutaric acid had a large decrease (− 2.15-fold) in eC. Ile was the only amino acid impacted in the eC plots, and was moderately increased by 0.81 fold (Table 1). Quinic acid, an important intermediate of the lignin pathway, increased in the eC plots (0.92-fold). 1-Benzylglucopyranoside had 0.95 log2-fold change in the eC plots.
The combination of eC and eT showed similar trends to the eT plots alone, but with a dampened effect (Table 1). All metabolites that had a significant change in the eCeT treatment were also significantly altered in the eT treatment. In both eT and eCeT, α-tocopherol increased significantly with a log2-fold > 1. Many compounds that were significantly increased in eT were not significantly altered in eCeT. For example, Ser and Thr were significantly increased in the eT treatment, but were not significantly impacted in eCeT, with a fold-change trend of decreasing in relative leaf content. Additionally, many sugars and sugar alcohols including; arabinose, fructose, galactose, glucose, and inositol increased in the eT plots but not in eCeT.

Integration of metabolite and transcript datatimepoint A
Weighted gene correlation network analysis (WGCNA) was used to construct a signed network of clustered genes from the time point A (sampled after 30 days of treatment). Using 9802 transcripts, expression data were clustered into 9 unique modules. GO enrichment was performed for each module to identify the overrepresented biological functions within each grouping (Supplemental Table 1). Module 1 contained the largest number of transcripts (1719 total), but only 2 were differentially expressed in the eCeT treatments. Module 3 contained 948 transcripts, 84 of which were differently expressed in eCeT, and had significantly enrichment in regulation of salicylic acid biosynthesis process, defense to oomycetes, cell surface receptor signaling, plant-type hypersensitive response, protein phosphorylation, defense response to bacterium, and oxidation-reduction process. In addition, module 3 had a significant negative correlation with melibiose (Supplemental Table 1, Supplemental  Fig. 4). Module 5 contained 484 transcripts, 35 of which were differently expressed in eCeT treatments, with significantly enrichment in oxylipin biosynthesis, terpenoid biosynthesis, response to wounding, hormone biosynthesis, secondary metabolite biosynthesis, defense response, and oxidation-reduction process. Within module 3, four transcripts were strongly correlated with melibiose content (r ≥ |0.8|, p < 0.01; Fig. 2). Melibiose showed significant negative correlations with the expression of two guard cell transcripts. The first, At4g23160, is annotated as cysteine-rich receptor-like protein kinase 8 (CRK8) and shows biological activity involved in defense response to bacterium and protein phosphorylation. The second was At3g48300, a cytochrome P450 (CYP71A23) involved in oxidation-reduction processes. The strongest relationship was identified for transcript ID, At1g54960, which was annotated as a subfamily of mitogen-activated protein kinase kinase kinase 2 (MAP3 K) known as Arabidopsis Nucleusand Phragmoplast-Localized Kinase1 (NPK1)-related protein kinases (ANPs) specifically ANP2, in the Arabidopsis genome. The final gene with a negative linear relationship with melibiose was identified as transcript ID, At5g54160, a caffeic acid/5-hydroxyferulic acid O-methyltransferase 1 (COMT) with functions in lignin and flavonoid biosynthesis.

Integration of metabolite and transcript datatimepoint B
Several significant relationships were identified within the metabolite and transcript datasets, with a r ≥ |0.75| and an adjusted p ≤ 0.01 (Table 2). All significant linear correlations between DEG and metabolites are shown in Supplemental File 3. Xylose had the most significant linear correlations, 10 positive and one negative. The strongest relationship was between gene ID At5g19440 (r = 0.87) (Fig. 3). The BlastX annotation hit for Arabidopsis thaliana identified At5g19440 as a NAD(P)-binding protein, located in the cytosol, and involved in cinnamyl-CoA reductase (CCR) activity in a variety of biological processes including lignin biosynthesis and response to cold. This sequence was highly similar to the annotated CCR1 in the Z. mays genome (LOC100382657) and predicated CCR1 in S.italica genome (LOC101786254). The second strongest relationship between a transcript and xylose was transcript At5g56870 (r = 0.84). This transcript was annotated as β-galactosidase 4 (BGAL4) in the Arabidopsis genome. This gene is localized to the cytoplasm and functions in rRNA N-glycosylase activity and defense response, also acting as a negative regulator of translation. Along with the strong linear correlation to xylose, BGAL4 was also significantly correlated with α-tocopherol, stigmasterol, and O-acetylsalicylic acid (r = 0.80, 0.76, 0.75, respectively; Supplemental File 3). Two additional transcripts, At5g03840 (r = 0.81) and At1g36160 (r = − 0.81), showed significant linear correlations with xylose content. Transcript ID At5g03840 was annotated as terminal flower 1 (TFL1) controlling inflorescence meristem identity  within the Arabidopsis genome and showed 73% homolog with the Zea mays ZCN12 (Zea centroradialis (ZCN)) (NM_001112779.2). This gene functions in many biological roles including acting as a negative regulator of cell aging and is expressed during growth and development stages from flowering to leaf senescence. ZCN12 was also found to have strong correlations with sinapic acid and valine (Supplemental File 3). Unlike the first three gene correlations, At1g36160, showed a negative linear correlation with xylose and sinapic acid (r = − 0.74). This gene was annotated in the Arabidopsis genome as acetyl-CoA carboxylase 1 (ACC1), which regulates the initial step in the fatty acid biosynthesis pathway of acetyl-CoA to malonyl-CoA. In addition to the strong associations with xylose (Fig. 3e), transcript ID At1g65480 was correlated with 4 additional metabolites. The BlastX annotation hit for Arabidopsis thaliana identified this transcript as flowering locus T (FLT) functioning in phosphatidylethanolamine binding with biological processes in flower development. The Blast annotation for Zea mays had only 52% homology with TWIN SISTER of FT (NM_001309850.1) and predicted Setaria italica FLT (XM_004960683.3) with a homology of 53%. Positive linear correlations were identified between FLT and valine, sinapic acid, O-acetylsalicylic acid, β-sitosterol, and xylose (r = 0.84, 0.82, 0.78, 0.77, 0.75) respectively (Supplemental Fig. 5).

Discussion
This study investigated the transcriptional and metabolite response of Guinea grass to growth at eC, eT and the combination of eCeT under open-air field conditions. Previous physiological analysis of the plants grown at this FACE experiment revealed that growth at eT increased leaf biomass and leaf area, but there was considerable variation within each treatment (Supplemental Fig. 6) (Britto de Assis Prado et al. 2016). The transcript and metabolite data also showed significant variation within each treatment, revealed by MDS (Fig. 1), but different groupings were apparent in the transcript and metabolite data. The transcriptome showed a clear separation among sampling dates (30 or 50 days of treatment), while the metabolome data showed a clear separation of treatments, most clearly C and eT (Fig. 1). The transcriptome also showed separation of treatments at the first sampling date (Timepoint A), but not the second (Timepoint B) (Fig. 1a). The withintreatment variability observed in the data may be due to environmental and edaphic properties of the field site (i.e., soil characteristics and/or slope), which could not be accounted for in the analysis, or from the relatively low sample size (n = 4). However, the fact that no transcripts were significantly altered in response to eC alone is consistent with previous experiments done with other C 4 grasses, which showed that in the absence of drought stress, there is no direct effect of elevated [CO 2 ] on C 4 metabolism (Wall et al. 2001;Leakey et al. 2006). On the other hand, eT stimulated the growth of P.maximum in this study (Supplemental Fig. 6), yet only a few genes were significantly differentially expressed in eT alone. Optimal growth for P. maximum is estimated to be at least 21 °C, and background temperature in ambient conditions during this experiment was only 17.0 ± 0.06 °C, so the eT treatment moved the plants closer to optimal growth temperatures (Araujo et al. 2013). In the combination (eCeT) treatment, significant changes in the transcriptome were observed, most notably transcripts involved in protein translation and assembly.
A key aim of this work was to integrate the transcript and metabolite data sets. Because of the clear differences in the transcript data between time points A and B (Fig. 1a), integration was done within each time point and not across time points. In addition, time point B offered the potential to integrate leaf physiological data with the transcript and metabolite datasets. Because of the large level of variability within individual treatments for all of the parameters, linear correlation analysis was used for data integration. For time point A WGCNA networks were built to group genes responding in similar ways together to form a larger network. These networks were then correlated with the metabolite data to identify interactions of interest. Because there were fewer DEG within time point B, WGCNA was not possible, and so for time point B, simple linear correlations between the DEG, significant metabolites, and leaf physiological data were used to identify correlations of interest.

Integration of metabolomic and transcript data-time point A
Two modules (3 and 5) from the WGCNA analysis contained significant numbers of DEG (Supplemental Table 2). The DEG identified within module 3 were enriched in salicylic acid (SA) biosynthesis (p 0.016) and general defense response to bacteria, suggesting greater stress in the eCeT treatment compared to control. The enrichment of SA biosynthesis transcripts has previously been identified to be part of the heat-stress response in plants. Salicylic acid stabilizes heat shock transcription factors increasing their binding to heat shock elements and thus promoting heat shock-related genes (Wahid et al. 2007). While the background temperature during the experiments was low, it is possible that the combination of eC and eT increased leaf temperatures more than just the eT temperature alone, thus triggering some stress-related transcriptional responses. The GO-term enrichment analysis of the transcripts associated with module 5 are also in agreement with a trend towards up-regulation of secondary metabolites involved in defense response. Many amino acids are precursors to secondary metabolites that change in response to various environmental perturbations. In this study, the eT and eCeT treatments resulted in significantly greater levels of amino acids derived from oxaloacetate and pyruvate. Accumulation of such amino acids is well known in abiotic stresses and heat (Joshi et al. 2010;Kaplan et al. 2004;Rizhsky et al. 2004;Guy et al. 2008;Obata and Fernie 2012;Obata et al. 2015). The role of branched chain amino acids and other amino acids in improving heat tolerance is complex and not fully understood, but likely provides support for an increase in production of secondary metabolites to influence defense response during temperature stress (Kaplan et al. 2004;Du et al. 2011). Additionally, increases in amino acid pools have the potential to serve as osmolytes under heat stress (Joshi et al. 2010) or an alternative electron donor for the mitochondrial electron transport chain (Araujo et al. 2013). High levels of Ile, Val, and Phe were observed in the leaves of Bermuda grass and mature kernels of rice exposed to elevated temperature and were thought to improve heat tolerance (Du et al. 51 Page 10 of 13 2011; Yamakawa and Hakata 2010;Kaplan et al. 2004). The buildup of transcripts in the brown and green modules involved in defense and secondary metabolite biosynthesis in combination with the increase in amino acids suggest upregulation of defense metabolism in the eCeT treatment at 30 days post-treatment exposure.
Plant grown under abiotic stress condition use compatible osmolytes, low molecular mass compounds, to provide an adaptive mechanism. Accumulation of these soluble sugars such as fructose, sucrose, and glucose are documented to signal stress conditions and improve tolerance to cold, heat, and drought stress ( (Du et al. 2011;Wahid et al. 2007;Kaplan et al. 2004;Rizhsky et al. 2004). The eT plots, regardless of [CO 2 ], showed increased soluble sugar content; including melibiose. Melibiose is a common osmolyte and a byproduct of the hydrolysis of raffinose, via invertase (ElSayed et al. 2014), and has been identified to accumulate under various environmental conditions including drought stress and elevated temperature in both C 3 and C 4 species (Rizhsky et al. 2004;Kaplan et al. 2004;Du et al. 2011). Figure 2 shows the relationship of four transcripts negatively associated with melibiose, the strongest of which was MAP3 K. MAPKs are highly conserved signaling modules in eukaryotes involved in cellular responses, including response to stimuli. These cascades trigger activation of transcription factors, via phosphorylation, and other cellular responses to signal a stimulus response (Xu and Zhang 2015). This transcript was annotated as ANP2/NPK1; these enzymes activate stress-related MAPKs activated by ROS (Kovtun et al., 2000). In this study, there was an inverse relationship between MAP3 K expression and melibiose content, suggesting that there was less stress-related signaling in the eT plots.
The second transcript highly associated with melibiose was a caffeic acid/5-hydroxyferulic acid O-methyltransferase 1 (COMT). Switchgrass and maize mutants have shown that downregulation of COMT impacts lignin content of biofuel crops specifically decreasing the S-lignin units and improving ethanol yields Guillaumie et al. 2008). Downregulation of this transcript could signal a shift from lignin biosynthesis to secondary defense metabolism.
The final two transcripts with a strong negative relationship to melibiose were CRK6 and CYP71A23, annotated as CYP72A26 in the maize genome (Jameson et al. 2008). Both of these transcripts are expressed in the guard cells. CRKs are induced by pathogen infection, and transgenic Arabidopsis plants with significantly elevated levels of CRK5 had enhanced leaf growth and increased pathogen resistance (Chen et al. 2003). CRK6 was found to be involved in ROS mediated signal processes with increased transcript abundance after exposure to pathogens, SA, and ROS (Chen et al. 2003(Chen et al. , 2004. The decrease in abundance of transcripts involved in guard cell metabolism is interesting and could potentially indicate a stomatal response within the eT and eTeC plots. The correlation between transcripts that are tightly regulated with complex feedback interactions between environment and cellular metabolism and melibiose lead us to suggest that melibiose functions in osmotic adjustment under eT and eTeC to protect cellular components and enhance acclimation to eT.

Integration of growth, metabolomic, and transcript data-timepoint B
Time Point B (50 d of treatment) allowed us to test the correlation of leaf biomass and area with metabolite and transcript levels; however, no significant correlations were observed. The low temperatures reported on the day of sampling may have greatly affected the metabolite and transcript levels, hindering the correlation of these datasets with leaf measurements taken at the end of the season. Field conditions are also inherently more variable than controlled environment conditions, but field experiments better represent the variation expected in the future. After 50 d of treatment, we did observe significant correlations between transcripts and metabolites, with the most notable being xylose (Table 2). Xylose is a monosaccharide and a precursor to pyruvate and acetyl-CoA. The use of linear correlation analysis allowed for the identification of associations that might not have been predicted with traditional methods. Three transcripts, shown in Fig. 3, had a strong positive correlation with xylose content. The strongest correlation was with Z. mays cinnamoyl-CoA reductase (CCR1). CCR1 is the entry point into the lignin branch of the phenylpropanoid pathway (Piquemal et al. 1998). The position of CCR in the phenylpropanoid pathway gives this enzyme direct control over the direction of the metabolic flux toward flavonoids or monolignols (Sattler et al. 2017). ZmCCR1 mutants with a moderate down-regulation of CCR1 activity increased the digestibility of cell walls without severely modifying lignin content (Tamasloukht et al. 2011). The overexpression of SbMyb60 transcription factor in Sorghum bicolor was found to directly impact the monolignol biosynthesis with higher abundance of syringyl lignin, cell wall composition, with increases in soluble phenolic and aromatic amino acids (Scully et al. 2016(Scully et al. , 2018. The observed increase in abundance of sinapic acid, phenylalanine, and α-Tocopherol under the eT treatment, regardless of [CO 2 ], point to an impact of eT on lignin composition, possibly increased S-type concentrations, although that would need to be tested. In addition to components of cell wall structure, monolignol polymers have been induced under various biotic and abiotic stress conditions (Tronchet et al. 2010;Vincent et al. 2005). A similar heating experiment with both C 4 and C 3 grasses was performed at the prairie heating and CO 2 enrichment site located in Wyoming, USA. A metabolomic profiling of Bouteloua gracilis, a C 4 grass, identified similar responses including an increase in phenolic and sugar content within the senesced leaves of plants grown under elevated temperature and the combination treatment, along with an increase in amino acid content under the elevated temperature with ambient [CO 2 ] (Suseela et al. 2014).
Sugars are continuously moving thru the phloem from source to sink and convey signaling information along the way, including regulation of flowering time (Gibson 2005). In several species adding exogenous sucrose and possibly other sugars promotes flowering (Cho et al. 2018). Interestingly, xylose was identified to have two significant positive correlations with flowering time-related genes, FT and ZCN12 (AtTFL1, Fig. 3c, d). The Arabidopsis orthologs are two closely related genes that are key integrators of flowering transition pathways. These two genes have antagonistic functions in the flower transition pathway; FT promotes flowering and TFL1 suppresses the transition (Shannon and Meekswagner 1991). Unlike the ortholog gene, ZCN12 belongs to the FT-like II group, primarily expressed in the leaf blades, and is activated after floral transition and is continually express thru late stages of vegetative development (Danilevskaya et al. 2010). FT-like proteins mediate numerous processes including growth, plant architecture, fruiting and tuber formation (Pin and Nilsson 2012). Recent evidence has shown FT can act as a cell autonomous regulator of stomatal guard cell opening and closure (Kinoshita et al. 2011). Kinoshita and colleagues showed overexpression of FT in Arabidopsis increased H + -ATPases activity within the guard cells, suggesting that FT likely influences the stomatal opening/closure via regulation of H + -ATPase. The accumulation of FT related transcripts within eT and eCeT could offer an additional regulation of stomatal function correlating with the osmotic adjustment provided from the accumulation of soluble sugars, such as melibiose and xylose.

Conclusion
This study investigated tropical Guinea grass response to global atmospheric change and rising temperature in the field. Although warming treatments in temperate regions often decreased productivity, heating the canopy 1.5 °C above ambient in this experiment resulted in increased leaf area and biomass. Field transcriptomics and metabolomics identified metabolic pathways that were altered by growth at elevated [CO 2 ] and temperature. Melibiose and xylose content were greater in plants exposed to elevated temperature, and content of these sugars was significantly correlated to the expression of genes involved in secondary metabolism, defense response, and stomatal function. We hypothesize that the metabolite and transcript responses were associated with alleviation of stress and greater growth at elevated temperature. 1238030 to EAA. We thank Amy Marshall-Colon for assistance with statistical analysis. Any opinions, findings, conclusions, or recommendations expressed in this publication are those of the author(s) and do not necessarily reflect the view of the U.S. Department of Agriculture. Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the U.S. Department of Agriculture. USDA is an equal opportunity provider and employer.
Author contributions CAM and EAA conceived and designed the study. JMW, CRY, TRM and SC performed the experiments. TRM and SC collected leaf samples from the field and did the RNA extraction. CRY did the RNAseq experiments and analyzed transcriptome data. JMW did the GC/MS experiments and analyzed the transcriptional and metabolite data and wrote the paper with EAA. All authors provided critical feedback on the research, analysis and manuscript.

Data availability
The metabolomics and metadata reported in this paper are available via the NIH Common Fund's Data Repository and Coordinating Center (supported by NIH grant, U01-DK097430) website, the Metabolomics Workbench, http://www.metabolomicsworkbench.org, where it has been assigned Project ID (PR000717). The data can be accessed directly via it's Project DOI: 10.21228/M80M50. The transcriptomic and metadata reported in this publication have been deposited in NCBI's Gene Expression Omnibus and are accessible through GEO series accession number GSE122194 (https ://www.ncbi. nlm.nih.gov/geo/query /acc.cgi?acc=GSE12 2194https://www.ncbi.nlm. nih.gov/geo/query/acc.cgi?acc=GSE122194).

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.