Integrated Transcriptomic and Metabolomics Analysis of the Root Responses of Orchardgrass to Submergence Stress

Submergence stress can severely affect plant growth. Orchardgrass (Dactylis glomerata L.) is an important forage grass, and the molecular mechanisms of orchardgrass to submergence stress are not well understood. The roots of the flood-tolerant cultivar “Dian Bei” were harvested at 0 h, 8 h and 24 h of submergence stress. The combined transcriptomic and metabolomic analyses showed that β-alanine metabolism, flavonoid biosynthesis, and biosynthesis of amino acid pathways were significantly enriched at 8 h and 24 h of submergence stress and were more pronounced at 24 h. Most of the flavonoid biosynthesis-related genes were down-regulated for the synthesis of metabolites such as naringenin, apigenin, naringin, neohesperidin, naringenin chalcone, and liquiritigenin in response to submergence stress. Metabolites such as phenylalanine, tyrosine, and tryptophan were up-regulated under stress. The predominant response of flavonoid and amino acids biosynthesis to submergence stress suggests an important role of these pathways in the submergence tolerance of orchardgrass.


Introduction
Submergence is a major environmental stress limiting plant growth and development. Plants subjected to complete or partial submergence stress often suffer from oxygen deficiency. It has been found that submergence stress affects plant growth and physiology, including but not limited to reducing shoot and root growth, leaf chlorophyll, carbohydrate content, photosynthetic rate, stomatal conductance, and antioxidant enzyme activities, as well as activating anaerobic fermentation and altering hormone interactions in different plant species [1][2][3][4][5][6][7][8]. The tolerant species or genotypes showed higher root dry weight, chlorophyll content, antioxidant capacity, soluble sugar content and protein levels, and a lesser degree of lipid peroxidation and lower electrolyte leakage in roots or shoots under submergence than the sensitive ones [9,10]. Previous research has indicated some crucial physiological responses associated with plant survival, but the adaptative mechanisms to submergence stress are not fully understood due to genotypic variations in response to stress and recovery from the stress.

Root Appearance after Submergence Stress
The root tips became darkened after exposure to submergence stress and were more pronounced at 24 h of stress ( Figure 1). In addition, adventitious roots were produced at the base of the stem at 24 h of stress, but root length was not significantly altered. Meanwhile, the above-ground parts tended to elongate as the duration of submergence increased. This suggested that short-term submergence stress already induced leaf elongation of the orchardgrass, which might help resist the extreme hypoxic environment.
integrating transcriptomic and metabolomic analyses. Therefore, we designed this e iment to explore transcriptomic and metabolomics alterations in the roots of orchard under submergence stress. The correlations between gene expression and metab were analyzed to gain a better understanding of the molecular responses of roots t stress. The results will be helpful for genetic improvements in orchardgrass aim breeding new varieties with enhanced submergence tolerance.

Root Appearance after Submergence Stress
The root tips became darkened after exposure to submergence stress and were pronounced at 24 h of stress ( Figure 1). In addition, adventitious roots were produc the base of the stem at 24 h of stress, but root length was not significantly altered. M while, the above-ground parts tended to elongate as the duration of submergen creased. This suggested that short-term submergence stress already induced leaf el tion of the orchardgrass, which might help resist the extreme hypoxic environment.

Transcriptome Sequencing Analysis and Identification of DEGs
The raw data were filtered, checked for sequencing error rates and GC conten tribution, and then clean reads (Table S1) were obtained for subsequent analysis. The reads were mapped to the orchardgrass reference genome using HISAT2 software ( S2). The results showed that the proportion of sequenced reads successfully mapp the genome was higher than 70%, indicating the high reliability of the sequencing The correlation of gene expression levels between samples showed good intra-grou producibility ( Figure S1).
Based on the sequencing data, we analyzed the RNA-Seq at 0 h, 8 h, and 24 h to r the DEGs under submergence stress. The distribution of DEGs in different treat groups is shown in Figure 2a

Transcriptome Sequencing Analysis and Identification of DEGs
The raw data were filtered, checked for sequencing error rates and GC content distribution, and then clean reads (Table S1) were obtained for subsequent analysis. The clean reads were mapped to the orchardgrass reference genome using HISAT2 software (Table S2). The results showed that the proportion of sequenced reads successfully mapped to the genome was higher than 70%, indicating the high reliability of the sequencing data. The correlation of gene expression levels between samples showed good intra-group reproducibility ( Figure S1).
Based on the sequencing data, we analyzed the RNA-Seq at 0 h, 8 h, and 24 h to reveal the DEGs under submergence stress. The distribution of DEGs in different treatment groups is shown in Figure 2a. A total of 6663 DEGs (3096 up-regulated and 3567 down-regulated) were identified in the roots after 8 h of submergence, while 9857 DEGs (4779 up-regulated and 5078 down-regulated) were identified at 24 h of treatment (Figure 2b,c). Among them, 5455 DEGs were co-expressed. The gene expression pattern changed significantly after submergence stress (Figure 2d).

GO Analysis and KEGG Pathway Analysis
To fully understand the functional properties of DEGs and the metabolic pathways of gene products, we performed functional enrichment analyses. Each group of DEGs was annotated into three categories: biological process; molecular function; and cellular component ( Figure S2). The Gene Ontology (GO) enrichment analysis revealed that most of the terms had more down-regulated DEGs than up-regulated DEGs and more DEGs at 24 h than at 8 h. The DEGs were highly expressed in terms of oxidative stress, coenzyme-binding oxidoreductase activity, acting on peroxide as an acceptor, antioxidant activity, etc., under stress. Unlike at 8 h, the genes related to peroxidase activity were more active at 24 h. Moreover, some highly expressed genes were involved in multiple regulatory pathways simultaneously under 8 h and 24 h of submergence stress (Figure 3a,b, Table S3).

GO Analysis and KEGG Pathway Analysis
To fully understand the functional properties of DEGs and the metabolic path of gene products, we performed functional enrichment analyses. Each group of DEG annotated into three categories: biological process; molecular function; and cellular ponent ( Figure S2). The Gene Ontology (GO) enrichment analysis revealed that m the terms had more down-regulated DEGs than up-regulated DEGs and more DEG h than at 8 h. The DEGs were highly expressed in terms of oxidative stress, coenz binding oxidoreductase activity, acting on peroxide as an acceptor, antioxidant ac etc., under stress. Unlike at 8 h, the genes related to peroxidase activity were more at 24 h. Moreover, some highly expressed genes were involved in multiple regu pathways simultaneously under 8 h and 24 h of submergence stress (Figure 3a  The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment was analyzed to show the effect of submergence stress on the metabolic pathways of DEGs. The top ranked 20 metabolic pathways were significantly enriched at 8 h and 24 h, compared to 0 h, including phenylpropanoid biosynthesis, plant-hormone signal transduction, the biosynthesis of amino acids, starch and sucrose metabolism, plant-pathogen interaction, and MAPK signaling pathway-plant (Figure 3c,d). The number of DEGs in each metabolic pathway was higher at 24 h than at 8 h (Table S4). The antioxidant-related metabolic pathways, such as ascorbate and aldarate metabolism and glyoxylate and dicarboxylate metabolism, were significantly enriched at 24 h of submergence.

qRT-PCR Validation
To verify the reliability of the transcriptome sequencing data, we performed qRT-PCR validation on the 11 screened DEGs. Correlation analysis of qPCR and RNA-Seq results for these 11 genes revealed that R 2 was 0.931 at 8 h and 0.806 at 24 h of submergence stress. The validation results were basically consistent with the expression level of RNA-seq sequencing results ( Figure S3). This indicated that the data obtained by transcriptome sequencing were reliable.
pathway was higher at 24 h than at 8 h (Table S4). The antioxidant-related metabolic pathways, such as ascorbate and aldarate metabolism and glyoxylate and dicarboxylate metabolism, were significantly enriched at 24 h of submergence.

qRT-PCR Validation
To verify the reliability of the transcriptome sequencing data, we performed qRT-PCR validation on the 11 screened DEGs. Correlation analysis of qPCR and RNA-Seq results for these 11 genes revealed that R 2 was 0.931 at 8 h and 0.806 at 24 h of submergence stress. The validation results were basically consistent with the expression level of RNAseq sequencing results ( Figure S3). This indicated that the data obtained by transcriptome sequencing were reliable.

Quality Control of Metabolomic Data
A total of 439 metabolites were obtained from a control group of root samples extracted from 0 h of submergence stress and 8 h and 24 h as treatment groups for untargeted metabolomic analyses. A hierarchical clustering heat map was constructed and annotated with Human Metabolome Database (HMDB) classification. The phenylpropanoidand polyketide-related metabolites had the highest number under submergence stress (Figures 4a and S4). Correlation analysis of the QC samples showed R 2 values all close to 1 ( Figure S5), indicating that the entire analysis process was stable and reproducible. Subsequently, principal component analysis (PCA) was performed, and there was a significant separation between 0 h, 8 h, and 24 h (Figure 4b). PC1 and PC2 explained more than 35.0% of the variability, and metabolites were mainly distinguished by PC1 (Figure 4b). Through partial least-squares discrimination (PLS-DA) analysis, treatment groups were separate from each other, with both R 2 and Q 2 close to 1 ( Figure S6), indicating a stable and reliable relationship between metabolite expression and sample categories.
with Human Metabolome Database (HMDB) classification. The phenylpropano polyketide-related metabolites had the highest number under submergence stre ures 4a and S4). Correlation analysis of the QC samples showed R 2 values all clo ( Figure S5), indicating that the entire analysis process was stable and reproducibl sequently, principal component analysis (PCA) was performed, and there was a cant separation between 0 h, 8 h, and 24 h (Figure 4b). PC1 and PC2 explained mo 35.0 % of the variability, and metabolites were mainly distinguished by PC1 (Figu Through partial least-squares discrimination (PLS-DA) analysis, treatment group separate from each other, with both R 2 and Q 2 close to 1 ( Figure S6), indicating a and reliable relationship between metabolite expression and sample categories.

Differentially Expressed Metabolites (DEMs) and Functional Enrichment
A total of 120 DEMs were obtained at 8 h of submergence compared to 0 h 5a). Of them, 64 were up-regulated and 56 were down-regulated ( Figure 5b   The DEMs were subjected to KEGG-enrichment analysis, and the top 20 enriched pathways were identified (Table S5). The top five enriched metabolic pathways after 8 h of stress were pantothenate and CoA biosynthesis, tryptophan metabolism, β-alanine metabolism, flavonoid biosynthesis, and phenylalanine, tyrosine and tryptophan biosynthesis; meanwhile, the zeatin biosynthesis pathway was also significantly enriched at 8 h ( Figure 6a). The top five metabolic pathways enriched after 24 h of treatment were propanoate metabolism, alpha-linolenic acid metabolism, ascorbate and aldehyde metabolism, flavonoid biosynthesis, and vitamin B6 metabolism (Figure 6b). Flavonoid biosynthesis was significantly enriched in the roots at both 8 h and 24 h of stress.

Differentially Expressed Metabolites (DEMs) and Functional Enrichment
Subsequently, the DEMs with p-value < 0.05 and ranked in the top 20 (from smallest to largest p-values) were screened for correlation analyses and categorical annotation to verify synergistic or mutually exclusive relationships between different metabolites ( Figure S7). Flavonoid biosynthesis-related metabolites naringenin chalcone and naringenin were significantly enriched at 8 h of submergence stress, while naringenin chalcone, naringenin, apigenin, and scutellarin in this class were enriched at 24 h of treatment (Table S6). All the compounds in this class showed extremely strong synergistic effects. Moreover, amino acid biosynthesis-related metabolites, such as tryptophan, gammaaminobutyric acid, and n-epsilon-acetyllysine, were also significantly enriched and showed synergistic relationships with each other (Table S6). However, at 8 and 24 h of stress, the metabolite jasmonic acid had mutually exclusive relationships with most of the metabolites except for dihydrojasmone and n-[(−)-Jasmonoyl]-(l)-isoleucine, which showed strong synergistic relationships with each other (Table S6) The DEMs were subjected to KEGG-enrichment analysis, and the top 20 enriched pathways were identified (Table S5). The top five enriched metabolic pathways after 8 h of stress were pantothenate and CoA biosynthesis, tryptophan metabolism, β-alanine metabolism, flavonoid biosynthesis, and phenylalanine, tyrosine and tryptophan biosynthesis; meanwhile, the zeatin biosynthesis pathway was also significantly enriched at 8 h ( Figure 6a). The top five metabolic pathways enriched after 24 h of treatment were propanoate metabolism, alpha-linolenic acid metabolism, ascorbate and aldehyde metabolism, flavonoid biosynthesis, and vitamin B6 metabolism (Figure 6b). Flavonoid biosynthesis was significantly enriched in the roots at both 8 h and 24 h of stress.
Subsequently, the DEMs with p-value < 0.05 and ranked in the top 20 (from smallest to largest p-values) were screened for correlation analyses and categorical annotation to verify synergistic or mutually exclusive relationships between different metabolites (Figure S7). Flavonoid biosynthesis-related metabolites naringenin chalcone and naringenin were significantly enriched at 8 h of submergence stress, while naringenin chalcone, naringenin, apigenin, and scutellarin in this class were enriched at 24 h of treatment (Table  S6). All the compounds in this class showed extremely strong synergistic effects. Moreover, amino acid biosynthesis-related metabolites, such as tryptophan, gamma-aminobutyric acid, and n-epsilon-acetyllysine, were also significantly enriched and showed synergistic relationships with each other (Table S6). However, at 8 and 24 h of stress, the metabolite jasmonic acid had mutually exclusive relationships with most of the metabolites

C orrelation and Enrichment Analysis between DEGs and DEMs
To further analyze the role of genes and metabolites in the regulation of orchardgrass under submergence stress, a nine-quadrant diagram was drawn to illustrate the correlation between genes and metabolites. Among them, metabolite abundance located in quadrants one, two, and four was higher than gene abundance, indicating that metabolites were up-regulated while genes remained unchanged or down-regulated. Metabolite abundance located in the sixth, eighth, and ninth quadrants was lower than gene abundance, indicating that genes were up-regulated but metabolites remained unchanged or down-regulated. Only genes and metabolites located in the third and seventh quadrants showed consistent differential expression patterns between DEGs and DEMs. A total of 234 metabolites were positively regulated by 2389 genes at 8 h of submergence, and a total of 286 metabolites were positively regulated by 2389 genes at 24 h of stress (Figure 7a,b). These results suggest that these changes in metabolite accumulation might be directly or indirectly regulated by the corresponding genes.  (Table S6).

C orrelation and Enrichment Analysis between DEGs and DEMs
To further analyze the role of genes and metabolites in the regulation of orchardgrass under submergence stress, a nine-quadrant diagram was drawn to illustrate the correlation between genes and metabolites. Among them, metabolite abundance located in quadrants one, two, and four was higher than gene abundance, indicating that metabolites were up-regulated while genes remained unchanged or down-regulated. Metabolite abundance located in the sixth, eighth, and ninth quadrants was lower than gene abundance, indicating that genes were up-regulated but metabolites remained unchanged or down-regulated. Only genes and metabolites located in the third and seventh quadrants showed consistent differential expression patterns between DEGs and DEMs. A total of 234 metabolites were positively regulated by 2389 genes at 8 h of submergence, and a total of 286 metabolites were positively regulated by 2389 genes at 24 h of stress (Figure 7a,b). These results suggest that these changes in metabolite accumulation might be directly or indirectly regulated by the corresponding genes.
The DEGs and DEMs were simultaneously subjected to KEGG pathway analysis. The results showed that 16 and 27 metabolic pathways were enriched at 8 h vs. 0 h and at 24 h vs. 0 h, respectively. Among them, two metabolites and 17 genes were significantly enriched to β-alanine metabolism, three metabolites and 30 genes to flavonoid biosynthesis, and two metabolites and 27 genes to phenylalanine, tyrosine, and tryptophan biosynthesis at 8 h of submergence stress (Figure 7c). Three metabolites and 114 genes were enriched for the biosynthesis of amino acids, two metabolites and 29 genes for the biosynthesis of flavonoids, two metabolites and 53 genes for the metabolism of cysteine and methionine, and two metabolites and 37 genes for the biosynthesis of phenylalanine, tyrosine, and tryptophan biosynthesis at 24 h of submergence (Figure 7d). In addition, the metabolic pathways of zeatin biosynthesis, pantothenate, and CoA biosynthesis were significantly enriched at 8 h and 24 h of stress.

Biosynthesis of Flavonoid and Amino Acids
Because the flavonoid biosynthesis and biosynthesis of amino acids pathways were significantly enriched through the integrated DEGs and DEMs analysis under submergence, the expression patterns and network interactions of the DEGs and DEMs annotated in these two pathways were analyzed. Most DEGs associated with flavonoid biosynthesis The DEGs and DEMs were simultaneously subjected to KEGG pathway analysis. The results showed that 16 and 27 metabolic pathways were enriched at 8 h vs. 0 h and at 24 h vs. 0 h, respectively. Among them, two metabolites and 17 genes were significantly enriched to β-alanine metabolism, three metabolites and 30 genes to flavonoid biosynthesis, and two metabolites and 27 genes to phenylalanine, tyrosine, and tryptophan biosynthesis at 8 h of submergence stress (Figure 7c). Three metabolites and 114 genes were enriched for the biosynthesis of amino acids, two metabolites and 29 genes for the biosynthesis of flavonoids, two metabolites and 53 genes for the metabolism of cysteine and methionine, and two metabolites and 37 genes for the biosynthesis of phenylalanine, tyrosine, and tryptophan biosynthesis at 24 h of submergence (Figure 7d). In addition, the metabolic pathways of zeatin biosynthesis, pantothenate, and CoA biosynthesis were significantly enriched at 8 h and 24 h of stress.

Biosynthesis of Flavonoid and Amino Acids
Because the flavonoid biosynthesis and biosynthesis of amino acids pathways were significantly enriched through the integrated DEGs and DEMs analysis under submergence, the expression patterns and network interactions of the DEGs and DEMs annotated in these two pathways were analyzed. Most DEGs associated with flavonoid biosynthesis were down-regulated in expression after submergence stress, while the metabolites naringenin, apigenin, naringin, neohesperidin, and naringenin chalcone, were up-regulated in expression, suggesting that orchardgrass may regulate the biosynthesis of flavonoidrelated substances by down-regulating related genes in response to submergence stress (Figure 8a1,a2). The genes involved in the biosynthesis of the amino acid pathway were more significantly up-regulated at 24 h of submergence, and the metabolites tryptophan and l-saccharopine were both up-regulated at 24 h of stress (Figure 8b1,b2). The network analysis of DEGs and DEMs annotated in these two metabolic pathways showed positive and negative regulations of genes and metabolites (Table S7a,

Discussion
Through the combined transcriptomic and metabolomics analyses, the post-transcriptional state of related gene expression can be explored to further understand the response mechanisms of orchardgrass to submergence stress. When the roots of orchardgrass were exposed to 8 h and 24 h submergence, genes involved in flavonoid biosynthesis and amino acid biosynthesis were significantly expressed, and the content of compounds, such as flavonoids and amino acids, also showed an increasing trend. The results indicate that the flavonoid and amino acid biosynthetic pathways largely responded to submergence stress in orchardgrass. Similar results were found in the aboveground tissue of red clover (Trifolium pratense) to Pb toxicity [34,35].
Flavonoids participate in plants' responses to environmental stress and serve as strong antioxidants and free radical scavengers [36]. A study on the radical scavenging and antioxidant activities of flavonoids from Scutellaria baicalensis extracts concluded that the scavenging activities of different flavonoids varied in response to different environmental conditions [37]. Submergence stress increased lipid peroxidation of roots and decreased antioxidant activities in four genotypes of perennial ryegrass under submergence stress [4]. In this study, naringenin, apigenin, naringin, neohesperidin, and naringenin chalcone involved in flavonoid biosynthesis were significantly enriched at 24 h of submergence stress. Flavonoid compound accumulation was also associated with the submergence tolerance of willow [25] and pot mum [22]. The results suggest the role of flavonoid in coping with oxidative injury from submergence stress.
It has been shown that the flavonoid biosynthetic pathway is regulated by different classes of transcription factors (TFs), such as MYB, WRKY, NAC and bHLH [38][39][40]. Among these TFs, the MYB-bHLH-WD40 ternary complex is one of the most important regulatory components involved in the biosynthesis of phenylpropanoids and flavonoids [22,41,42]. In this study, the metabolic pathway of flavonoid biosynthesis was enriched at 8 h and 24 h of submergence stress, accompanied by a high expression of genes encoding the TF family, such as WRKY and MYB. The results indicate that these TFs could regulate the biosynthesis of flavonoid compounds in orchardgrass roots under the water.
Amino acids are important physiologically active organic compounds in plants and play an important role in regulating growth. When plants experience abiotic stress, free amino acids can respond rapidly by increasing the concentration of certain amino acids for adaptation to environmental changes [43,44]. In this study, biosynthesis of amino acids and β-alanine metabolism were enriched at both the gene and metabolite levels in the roots of orchardgrass in response to submergence. Similar to the results of orchardgrass root under submergence stress, alpha-linolenic acid metabolism, phenylalanine metabolism, and phenylpropanoid biosynthesis pathways were also significantly enriched in the leaf transcriptome analysis [32]. Meanwhile, oxidoreductase activity, cofactor bling, and oxidationreduction processes are highly present in GO-enriched terms, suggesting that the leaves of orchardgrass have a complex oxidoreductase process in response to flooding stress. However, carbon metabolism, flavonoid biosynthesis, and zeatin biosynthesis pathways were not significantly enriched in the leaves [32]. The up-regulated expression of genes encoding β-alanine metabolism, methionine metabolism, tyrosine metabolism, glycine, serine, and threonine metabolism was also found in field cress (Sesbania cannabina) roots after 3 h of submergence stress [45]. In wheat, most of the measured amino acids increased in shoots during the first 12 days of submergence, with only five showing decreasing or unchanged levels, including alanine in both tolerant and sensitive cultivars. However, all amino acid levels were lower in the intolerant cultivar than in the tolerant cultivar on days 14 and 16 of treatment [24]. All results suggest that amino acid metabolism-related processes, such as β-alanine metabolism and phenylalanine metabolism, are associated with submergence tolerance, but the responses of amino acids to submergence varies with plant species, stress duration and intensity. We speculated that the increased amino acids in the roots of orchardgrass after submergence stress was probably due to accelerated nitrogen metabolism to make the raw material for more abundant protein biosynthesis.
The metabolic pathways of zeatin biosynthesis, pantothenate, and CoA biosynthesis were significantly enriched in orchardgrass at 8 h and 24 h of submergence stress, suggesting a role of these pathways in submergence tolerance. Zeatin is a natural cytokinin present in plants, which promotes plant-cell division, prevents chlorophyll and protein degradation, maintains cell viability, and delays plant aging [46,47]. The up-regulated zeatin biosynthesis found in this study indicated its role in promoting growth and assisting submergence tolerance in orchardgrass. The enriched DEGs for zeatin biosynthesis were also found in Populus ussuriensis exposed to salt stress [48]. CoA is an essential cofactor in the metabolism and is particularly important in the metabolism of fatty acids. A combined transcriptomic and metabolomic analysis of the response of Zygophyllum plants to salt stress revealed that pantothenate and CoA biosynthesis were significantly enriched in both salt-tolerant and sensitive varieties [49], demonstrating its role in altering plant response to stress conditions.

Plant Materials, Submergence Treatment and Harvesting
The flood-tolerant orchardgrass variety 'Dian Bei' was used for this experiment. The seeds were obtained from the College of Grassland Science and Technology, Sichuan Agricultural University, Chengdu, China. The seeds were placed in Petri dishes containing two layers of moist filter paper and germinated in dark conditions in an incubator at 22 • C. After germination, the uniformly growing plants were transplanted into pots (15.

Transcriptome Sequencing and Data Analysis
At 0 h, 8 h, and 24 h of submergence stress, approximately 5 g of the root system was taken for each sample, rapidly placed in liquid nitrogen, and then stored in a freezer at −80 • C. Three biological replicates were made at each time point. Total RNA was extracted using a Trizol kit (Invitrogen, Carlsbad, CA, USA), and the samples were tested for RNA integrity and purity by agarose gel electrophoresis using a Nano Photometer spectrophotometer (Implen, Westlake Village, CA, USA) and an Agilent 2100 bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The samples that met the requirements were sequenced by Illumina PE150 (Illumina, San Diego CA, USA) at Beijing NovoMagic.
The sequencing raw data were filtered by removing low-quality reads and detecting sequencing error rate and GC content distribution. The clean reads were mapped to the orchardgrass reference genome [50,51] using HISAT2 software (http://ccb.jhu.edu/ software/hisat2/faq.shtml, accessed on 20 July 2021) to obtain gene-on-gene localization information, as well as information on sequence characteristics specific to this species. Gene alignment statistics and fragments per kilobase of exon per million mapped fragments (FPKM) values were calculated. Differential expression results were analyzed for each comparison of treatment time points using DESeq2 (version 1.16.1) software [52], and padj < 0.05 and |log2FC| ≥ 1 were used as screening criteria for identifying the DEGs. To further implement the functions of DEGs, we used cluster profiler software to perform GO functional enrichment and KEGG-pathway-enrichment analysis.

qRT-PCR Analysis
Eleven DEGs with expression in both 8 h and 24 h of submergence stress were randomly selected for qRT-PCR using the orchardgrass actin gene as an internal reference gene [53,54]. Primers are listed in Table S8. RNA was reverse-transcribed into cDNA using a PrimeScript ™ RT reagent kit with a gDNA Eraser kit (Takara, Dalian, China), and qRT-PCR was performed with a TB Green Premix Ex TaqTM II kit (Takara, Dalian, China). Three technical replicates were performed for each target gene. The total volume for qRT-PCR was 10 µL with the following procedure: one cycle at 95 • C for 30 s; 40 cycles at 95 • C for 5 s and 60 • C for 34 s; then, 95 • C for 15 s; 60 • C for 1 min; and 95 • C for 15 s. Finally, the relative gene expression was calculated using 2 −∆∆Ct [55].

Metabolomics Data Analysis
Root tissues for metabolomics at different times of sampling were the same batch as those used for the transcriptome. A total of six biological replicates of each treatment were used for metabolomics analysis. Approximately 0.1 g of ground root tissue was placed in a tube with 500 µL of 80% aqueous methanol solution. The extraction was shaken and placed in an ice-water bath for 5 min, and then centrifuged at 15,000 r for 20 min at 4 • C. A certain amount of supernatant was taken and diluted with mass-spectrometry-grade water to make a solution containing 53% methanol. The supernatant was collected after centrifuging again for 20 min and used for liquid chromatography coupled with mass spectrometry (LC-MS) [56]. Briefly, an equal volume of samples was mixed as QC samples and 53% aqueous methanol solution was used as blank samples. The chromatographic conditions were as follows: column, HypesilGoldcolumn (C18), column temperature was 40 • C; flow rate was 0.2 mL/min; mobile phase A in positive ionization mode was 0.1% formic acid; mobile phase B was methanol; mobile phase A in negative ionization mode was 5 mM ammonium acetate (pH 9.0); and mobile phase B was methanol. The elution gradient is shown in Table S9. The mass spectrometry conditions were as follows: spray voltage of 3.5 kV; sheath gas flow rate at 35 psi; aux gas flow rate at 10 L/min; capillary temperature at 320 • C; S-lens RF level at 60; aux gas heater temperature at 350 • C; positive and negative polarity; and MS/MS secondary sweep at 350 • C. The MS/MS secondary scans were data-dependent scans. The reagents and instruments used in the experiment are shown in Tables S10 and S11.
The obtained raw data were imported into CD3.1 software, and each metabolite was screened for charge-to-mass ratio and retention-time parameters. The molecular formula of the metabolite was predicted by ion peaks and compared with the database. The relevant interfering ions were removed using blank samples, and the data were processed to obtain quantitative analysis results. Multivariate statistical methods, such as PCA and PLS-DA, were applied to reveal the differences in metabolic patterns among different treatments. Finally, DEMs were defined by VIP > 1.0, fold change (FC) > 1.2 or FC < 0.833, and p-value < 0.05. The biological significance associated with the metabolites was explained by functional analyses, such as metabolic pathways.

Integration of Transcriptomic and Metabolomic Analysis
The relationship between DEGs and DEMs were revealed based on Pearson's correlation coefficient analysis. All the obtained DEGs and DEMs were mapped to the KEGG pathway database to identify biochemical pathways and signal transduction pathways in which DEGs and DEMs were jointly involved. The common highly expressed pathways were identified for further analysis of the expression patterns and network interactions.

Conclusions
The results elucidated the molecular mechanisms of submergence tolerance in orchardgrass by exploring expressions of genes and metabolites in the roots under short-term stress conditions. The integrated transcriptomic and metabolomic analysis identified key metabolic pathways in response to stress. The up-regulated flavonoid biosynthesis and the up-regulated amino acids biosynthesis observed in the roots could contribute to submergence tolerance. These pathways could be potential targets for further explorations of the role of these pathways in submergence tolerance. Given the complex nature of submergence stress and potential genotype by environment interaction, the results would be helpful in designing future experiments for the dissection of these pathways in more diverse lines in response to different levels of submergence stress.