Metabolomic analysis with GC-MS to reveal potential metabolites and biological pathways involved in Pb & Cd stress response of radish roots

The radish (Raphanus sativus L.) is an important root vegetable crop. In this study, the metabolite profiling analysis of radish roots exposed to lead (Pb) and cadmium (Cd) stresses has been performed using gas chromatography-mass spectrometry (GC-MS). The score plots of principal component analysis (PCA) and partial least squares-discriminate analysis (PLS-DA) showed clear discrimination between control and Pb- or Cd-treated samples. The metabolic profiling indicated Pb or Cd stress could cause large metabolite alteration mainly on sugars, amino acids and organic acids. Furthermore, an integrated analysis of the effects of Pb or Cd stress was performed on the levels of metabolites and gene transcripts from our previous transcriptome work in radish roots. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of integration data demonstrated that exposure of radish to Pb stress resulted in profound biochemical changes including carbohydrate metabolism, energy metabolism and glutathione metabolism, while the treatment of Cd stress caused significant variations in energy production, amino acid metabolism and oxidative phosphorylation-related pathways. These results would facilitate further dissection of the mechanisms of heavy metal (HM) accumulation/tolerance in plants and the effective management of HM contamination in vegetable crops by genetic manipulation.

Heavy metal (HM)-pollution represents a major environmental hazard to human health due to the accelerating industrial developments in recent decades 1,2 . Among the non-nutrient heavy metals, lead (Pb) and cadmium (Cd) are the most common and widespread contaminants, which can be absorbed by plants through contaminated soil and water and then enter into the food chain 3,4 . In plants, heavy metal toxicity such as Pb and Cd could induce a number of morphological, physiological and biochemical defects including transpiration and photosynthesis alteration, carbohydrate metabolism imbalance and production of secondary stresses like nutrition stress, water deficit and oxidative stress [4][5][6] . Nevertheless, plants also have evolved diverse endogenous mechanisms to cope with the toxic effects of heavy metals including generating signal sensing and transduction proteins, activating transport systems and biosynthesis of chelating compounds [7][8][9] . However, systematically understanding the heavy metal toxicity mechanism and identifying the cellular or biochemical targets underlying plant physiological responses are complex and challenging.
Recently, the arrival of 'omics' approaches have been widely used in modern biology aimed at massively characterising the molecular mechanisms of living systems at different levels, which can give a system-wide view of understanding the various layers of the cellular architecture underlying function of biological networks 10-12 . Metabolomics is one of the last additions to the 'omics' wave, which can provide an approach to unravel the complex mechanisms by measuring many metabolites participating in various biochemical processes and across many biological systems 13,14 . There are various analytical platforms including gas chromatography-mass spectrometry (GC-MS), liquid chromatography (LC)-MS, capillary electrophoresis (CE)-MS and nuclear magnetic resonance

Results
Metabolic changes in response to Pb and Cd stresses in radish. In our preliminary experiment studies, different concentrations of Pb(NO 3 ) 2 (100, 200, 400, 800, 1,000 and 1,500 mg•L −1 ) and CdCl 2 •2.5H 2 O (50, 100, 150, 200, 400 and 500 mg · L −1 ) were set to investigate the changes of visible physiology symptoms in an advanced inbred radish line 'NAU-RG' with different temporal durations. Interestingly, no obvious morphologic differences were found among individuals when exposed to low dose of heavy metal treatment (0 to 800 mg · L −1 Pb(NO 3 ) 2 and 0 to 200 mg · L −1 CdCl 2 •2.5H 2 O) for a maximum of 72 h, while the plants were seriously hampered and grew abnormally when exposed to 1,500 mg · L −1 Pb(NO 3 ) 2 and 500 mg · L −1 CdCl 2 •2.5H 2 O. Therefore, we selected the concentrations of Pb(NO 3 ) 2 at 1,000 mg · L −1 (Pb1000) and CdCl 2 •2.5H 2 O at 400 mg · L −1 (Cd400) for the metabolomic analysis. Additionally, a control group was defined using non-treated seedlings (CK). Metabolites were extracted from taproot samples in eight replicate pools of plants for each of the three experimental groups including Pb1000, Cd400 and CK and were analyzed by GC-MS. The total ion chromatograms (TICs) of all samples (except one sample from Pb1000 group) demonstrated a strong signal, large peak capacity, and reproducible retention time, indicating the reliability of metabolomic analysis. Obvious chromatographic differences were observed between different sample groups, and a total of 1,104 types of metabolites were identified. Typical TICs from each group were shown in Fig. 1 and Supplementary Fig. S1. In order to reduce the dimensionality of the data and visualize samples grouping, an unsupervised multivariate data analysis method PCA was performed on the GC-MS data generated from Pb1000, Cd400 and CK groups. A PCA model was created with four principal components (PCs), which had explanation and predictability values of 74.8% and 51.9%, respectively ( Table 1). The score plot of the first two PCs was shown in Fig. 2A. Most of the data were lying inside the 95% confidence region (Hotelling T 2 ellipse). The results of the PCA for the three group samples indicated that an obvious separation between the control (CK) and treated samples (Pb1000 and Cd400) was detected, while no clear difference was observed between the Pb-(Pb1000) and Cd-treated (Cd400) groups ( Fig. 2A). In order to confirm this tendency, every two of these three groups were analyzed by PCA with similar results obtained. According to the PCA models, two, three and two PCs were gained from the control and Pb-treated, control and Cd-treated and Pb-and Cd-treated samples comparisons, respectively. The R 2 X and Q 2 values (goodness of prediction) are shown in Table 1, indicating that all the models could predictably explain the differences between groups. However, no distinct boundary between every two groups could be identified in the PCA score plots (Fig. 2B-D).
Compared with PCA, PLS-DA is a supervised method, which could classify the observations into the group from giving the largest predicted indicator variable. From that, two PCs (R 2 X = 0.560, R 2 Y = 0.516, Q 2 = 0.343) were indicated among the three groups, and the classification among these three groups was improved in the  Table 1. Explanation and predictability values of the principal component analysis (PCA) and partial least squares-discriminate analysis (PLS-DA). Control-Pb-Cd, analysed using the Pb-and Cd-treated samples as well as the controls; Control-Pb, analysed using the control and Pb-treated samples; Control-Cd, analysed using the control and Cd-treated samples; Pb-Cd, analysed using the Pb-and Cd-treated samples. score plot (Fig. 3A). More satisfactory modeling and prediction results with two PCs were obtained when data were analyzed only using the control and the Pb-or Cd-treated samples (R 2 Y > 0.9, Q 2 > 0.6) ( Table 1), and both Pb-and Cd-treated taproot samples were clearly separated from the control along the PC1 (Fig. 3B,C). As shown in Fig. 3D, Pb-and Cd-treated samples could be separated in the PLS-DA score plot with two PCs, despite minor overlap. Nevertheless, the R 2 Y and Q 2 values were only 0.476 and 0.092, respectively (Table 1), indicating a less metabolic change between these two treated samples when compared with the control.
Identification of the Pb-or Cd-responsive metabolites in radish roots. The altered metabolites were found from the line plots of the X-loadings of the first component of the PLS-DA pairwise comparison models. It was reported that the variable importance in the projection (VIP) values greater than 1 were considered the most relevant metabolites for explaining the responses 27 . On the basis of the parameter VIP > 1, 29 Pb-and 34 Cd-responsive metabolites with significant changes (student's T-test P < 0.05), were respectively identified (Supplementary Table S1 and Supplementary Table S2). Among them, three compounds were unknown, and one compound was a type of steroid that could not be accurately identified. The changed metabolites were mainly sugars, amino acids, and organic acids (Tables 2 and 3). Compared to the non-treated control, 10 metabolites increased and 18 decreased in levels of the radish roots after Pb exposure, while only one increased and 32 metabolites decreased by Cd treatment. Consistent with the results of PLS-DA score plots, a different but overlapping metabolite-response trend to Pb and Cd stress was determined in radish. Furthermore, the altered metabolites involved in the primary metabolism of glycolysis and citric acid cycle and the linking metabolites to amino acid synthesis in radish roots under Pb or Cd stress were respectively summarized on a simplified metabolic map (Figs 4 and 5). As shown in Figs 4 and 5, the accumulation of gluconate was both enhanced under the two treatments, while coordinated decreases in the content of five sugars (maltose, inositol, turanose, α -D-glucopyranoside and β -D-glucopyranose), two amino acids (alanine and glycine), three organic acids (acetimidic acid, decanedioic acid and oxalic acid), glycerol, hydroxylamine, phosphoric acid and sitosterol were identified. However, the contents of three sugars (fructose, galactose and glucose), three organic acids (citrate, hexadecanoic acid and octadecanoic acid), and a type of steride were increased in response to Pb treatment but decreased to Cd treatment. Furthermore, induced accumulation of α -linolenic acid and declined contents derived from malate, serine and isoleucine were identified by Pb treatment, while no significant difference in the contents of these metabolites was detected after exposure of Cd. In contrast, the levels of α -D-glucopyranoside, proline, threonine, pyroglutamate, 4-aminobutyrate, pyrophosphate, nicotinic acid, cholesterol, monostearin and two unknown compounds were decreased by Cd rather than Pb exposure. Pathway mapping and the metabolite-to-metabolite network visualization. All of the changed metabolites affected by Pb and Cd stresses were mapped to the biological pathways involved in the KEGG database, which were assigned to 49 pathways in either treatment (Supplementary Tables S3 and S4). The most statistically enriched pathways were analyzed with a Bonferroni correction (P ≤ 0.05). The results showed that nine and four pathways were enriched with changed metabolites, as a result of Pb and Cd exposure, respectively (Table 4). Among them, three, including galactose metabolism, starch and sucrose metabolism, and aminoacyl-tRNA biosynthesis, were enriched for both the Pb-and Cd-regulated metabolites. Furthermore, using all the altered metabolites as inputs, we constructed the metabolite-to-metabolite interaction networks comprising 14 and 16 metabolites for the Pb-and Cd-stress exposure in radish roots, respectively ( Supplementary Fig. S2). Both the networks could be divided into two sub-clusters. One sub-cluster consisted of compounds mainly involved in carbohydrate and glycan biosynthesis and metabolism, and the other was composed of several kinds of acidic compounds including amino acids, organic acids, and inorganic acids.

Integrated transcriptome and metabolome analysis of radish in response to Pb or Cd stress.
To gain a deeper understanding of the molecular mechanism of heavy metal tolerance and homeostasis in radish, an integrated analysis of the effects of Pb or Cd stress was performed on the levels of metabolites and gene transcripts from our previous transcriptome work in radish roots 24,25 . All the Pb-or Cd-responsive genes and metabolites were separated, and a gene-to-metabolite network was constructed each for the Pb- (Fig. 6) and Cd-stress exposure (Fig. 7). During Pb exposure, most of intersected pathways between differentially expressed genes and altered metabolites were involved in carbohydrate metabolism (i.e., starch and sucrose metabolism, fructose and mannose metabolism and the pentose phosphate pathway), energy metabolism (i.e., citrate cycle and glycolysis/gluconeogenesis), and glutathione metabolism (Supplementary Table S5). The gene-to-metabolite network consisted of 18 metabolites and 58 genes ( Fig. 6 and Supplementary Table S6). As shown in Fig. 6, the expression of several genes encoding enzymes including glutathione S-transferase (GSTU3, GSTU8, GSTU11 and GSTU23), carbonate dehydratase (ACA7 and BCA3), amidohydrolase (AAH), adenine phosphoribosyl transferase (APT3), glucose-6-phosphate dehydrogenase (G6PD2), phosphoserine aminotransferase (PSAT), adenosine deaminase (ADA), and nucleoside triphosphatase (APY1) were correlated with the regulation of glycine The fold changes in the concentrations of each metabolite between the two groups were calculated using the formula log 2 (Pb-treated/Control). "*" and "**" indicate the difference is significant (P < 0.05) and highly significant (P < 0.01) compared to the control, respectively. VIP, variable importance in the projection. content. Isopropylmalate isomerase-and tyrosine aminotransferase-encoding genes (IPMI2 and TyrAT) were found to be responsible for the decrease of isoleucine after the treatment of Pb. Six genes simultaneously participated in the regulation of alanine and serine, while three aminotransferase (HMT3, TyrAT and SERAT2;1), one cysteine synthase (OASC), one pyrimidine (PYD4), and one fructose-bisphosphate aldolase (FBA2) gene showed involvement. Additionally, most other genes including those encoding amylase (AMY1), 1,4-beta-D-xylan synthase (CSLD5), galacturonosyltransferase (GAUT14 and GAUT15), trehalose-phosphate synthase (TPS1 and TPS9), UDP-D-glucuronate 4-epimerase (GAE6), UDP-glucose 6-dehydrogenase (UDG4), galactosidase (AGAL1), fructose-bisphosphate aldolase (FBA2), phosphoenolpyruvate carboxykinase (PCK2), sucrose synthase (SUS5), glucan phosphorylase (PHS1), fructose-2,6-bisphosphatase (F2KP) and phosphofructokinase (PFK), Chitinase-like and G6PD2 were involved in the regulation of the concentration of storage carbohydrate (such as fructose, glucose, maltose and glucopyranose). All of the intersected pathways in response to Pb exposure were included in those pathways in response to Cd exposure. Also, some specific metabolites and genes responding to Cd stress were involved in pathways participating in oxidative phosphorylation, amino acid-related metabolism, and nitrogen metabolism (Supplementary Table S7). The corresponding gene-to-metabolite network under Cd exposure was mapped in Fig. 7, which included 23 metabolites and 82 genes (Supplementary Table S8). Similar genes were found to regulate the accumulation of inositol (i.e., NPCs, SAL2, PFK3 and ATP14K-α) and cholesterol (i.e., CYP707A1 and SMT2) in response to both Pb and Cd stress. Five GSTU genes, SHM7 (serine hydroxymethyltransferase-encoding gene), and 6PGD (6-phosphogluconate dehydrogenase-encoding gene) were reported to affect the accumulation of both glycine and pyroglutamate. Genes PYD4, PPOX, HEMA2, MTO2, SHM7, PSAT, TSB and AK-HSDH were responsible for the regulation of glycine and threonine. In addition, AAH, ACA7, SIN-like, SUVH3, UGLYAH, FAC1 and GDH1 were involved in glycine metabolism. Among multiple genes involved in the Cd-altered alanine metabolism, there  The fold changes in the concentrations of each metabolite between the two groups were calculated using the formula log 2 (Cd-treated/Control). "*" and "**" indicate the difference is significant (P < 0.05) and highly significant (P < 0.01) compared to the control, respectively. VIP, variable importance in the projection. were five genes encoding aminotransferases and five genes involved in the ethylene synthesis pathway. Similar to the situation in the Pb-related gene-to-metabolite network (Fig. 6), complex relationships were also observed in the Cd-regulated storage carbohydrate metabolism. In addition to PFK3, AGAL1, ASD1, TSP9, CSLD5, GAUT14, AMY1, FBA2, chitinase-like, PCK2 and MEE51, which also respond to Pb, genes encoding 6-dehydratase (RHM3), UDP-glucose (UDG4 and UGE3), GDP-L-fucose synthase (GER2), trehalose-6-phosphate phosphatase (TPPG), glucosidase (BGLU40), and ribose 5-phosphate isomerase (RPI2) were also responsible for Cd-regulated carbohydrate metabolism.

Discussion
Radish, an important vegetable crop, can uptake and accumulate heavy metals (HMs) including Pb and Cd in the root 28,29 . Many previous studies have reported the defense and damage of radish response to HMs at morphological, physiological and transcriptional levels [30][31][32][33] . However, there are no reports to investigate the metabolite changes of radish under the HMs exposure using metabolomics technology, which has been regarded as an important research field in the post-genomic area, especially for plant physiological responses to various stresses [34][35][36] . In this study, we report a comprehensive analysis of metabolic changes in radish responding to heavy metal (Pb or Cd) stress using a GC-MS-based metabolomics approach. The results indicated that heavy metal stress could cause many metabolite alterations, mainly alterations of sugars, amino acids, and organic acids. A KEGG pathway analysis integrating transcriptomic and metabolomic data demonstrated that exposure of radish to Pb stress resulted in profound biochemical changes including carbohydrate metabolism, energy metabolism and glutathione metabolism, while the treatment with Cd caused significant variation in energy production, amino acid metabolism and oxidative phosphorylation-related pathways. This study is the first report on the systematic identification and characterization of metabolic changes in response to HM stress in radish using metabolomic analysis.

Metabolic responses of radish to Pb and Cd stresses. Physiological processes such as photosynthesis
and respiration have been shown to be very sensitive to heavy metals in higher plants 37 . In our previous studies, chlorosis and growth inhibition has been observed as a primary toxicity symptom for radish under the high levels of heavy metals exposure such as Pb and Cd 30,33 , which may be attribute to the consequence of inhibited synthesis of chlorophyll, increased chlorophyll degradation and altered the levels of carbohydrate in the presence of metals. However, plants also have evolved intricate mechanisms to respond and defend the HM toxicity. For example, it was reported that the photoassimilates were stored in the form of hexoses and complex sugars under the HM stress, which can act as the osmoprotectants playing critical roles in osmotic adjustments or protection of cell

Pathway Adjusted P-value Count
Pb-treated stress constituents 21 . In the present study, the accumulations of glucose, galactose and fructose were enhanced under Pb treatment, while the contents of maltose, turanose, α -D-glucopyranoside and β -D-glucopyranose were reduced in the radish roots, suggesting that photoassimilates were stored as hexoses in radish after exposure to Pb. As shown in Fig. 4, the carbon flux mainly concentrated in the high diffusion rate of sugars (the upregulation of galactose, glucose and fructose) and glycolysis. The results may be attributed that the Pb stress has a negative effect on respiration 38 , therefore the radish roots need to strengthen the glycolysis to cope with the stress.
Scientific RepoRts | 5:18296 | DOI: 10.1038/srep18296 Figure 6. The gene-to-metabolite networks involved in radish root response to Pb stress. The annotations for the related Pb-responsive genes were listed in Table S6.  Table S8. Chelation of toxic heavy metal ions by low molecular weight organic acids (LMWOAs) is another commonly effective mechanism for the HM defense system 39 . Many studies have demonstrated that citrate, malate and oxalate could play central roles in case of tolerance to aluminum (Al) stress [40][41][42][43] . Additionally, citrate, oxalate and gluconate were observed to be significantly proficient under the Cd exposure 44 , and a positive correlation was also established between Pb and citrate concentrations in xylem sap of Sesuvium portulacastrum and Brassica juncea 45 . In this study, the altered organic acids involved in the primary metabolism of radish roots under Pb or Cd stress (Figs 4 and 5) were those including gluconate, citrate and malate. The content of gluconate was both enhanced under the two treatments, while up-regulated citrate and down-regulated malate were only shown during the Pb stress. In addition, some other organic acids which were not involved in the primary metabolism were significantly accumulated such as hexadecanoic acid, linoleic acid, α -linolenic acid and octadecanoic acid. These results indicated that LMWOAs were associated with HMs stress in case of radish. Particularly under the stimulus of Cd stress, carbon flow mainly concentrated in gluconate which revealed that the radish roots may be more inclined to use gluconate to alleviate the Cd toxicity (Fig. 5).
Many plants have been shown to accumulate proline and other amino acids when exposed to heavy metals, which could affect synthesis and activity of the antioxidant enzymes, regulation of ion transport, modulating stomatal opening and sequestering toxic cadmium ions in the cytosol or vacuole 21,46 . However, in our study, the content of all the altered amino acids were both decreased and the carbon fluxes were less flow to the process of amino acid synthesis upon the Pb (Fig. 4) or Cd exposure (Fig. 5) in radish. The result revealed that high level of Pb or Cd inhibited the synthesis of amino acids, which may be associated with the reduced activity of antioxidant enzymes (such as SOD, POD, CAT and APX) under the high concentrations of metal ions exposure reported from our previous studies 30,33 . Gene-to-metabolite networks of radish in response to Pb and Cd stresses. Integration of transcriptomics and metabolomics data can better elucidate the gene functions and metabolic pathways in biology 47 . In the current study, the gene-to-metabolite networks of radish responding to Pb and Cd stress were respectively constructed, which exhibited a complex regulation mechanism such as one compound was affected by several genes and one gene regulated many metabolites. In addition, the primary metabolic and transcript changes were summarized on a simplified metabolic map ( Supplementary Figs S3 and S4) according to the results of overlapping pathways involved in the KEGG metabolic database, which could strongly exhibit the direct metabolic relationship upon Pb and Cd exposure in radish roots. As shown in Supplementary Fig. S3 and S4, the expression of genes involved in glycolysis such as PFK3, FBA2 and PGK1 were affected by both the treatments of Pb and Cd stress. However, the contents of their direct products β -D-fructose-1, 6 P2, glyceraldehyde-3P and glycerate-3P were not influenced, suggesting that a complex regulation mechanism may exist.
Many studies reported that the glutathione (GSH) metabolism could play critical roles in heavy metal stress tolerance 48,49 . The GSTU gene families were the plant-specific tau class encoding glutathione S-transferases (GST), which were verified to participate in the conversion from glutathione and glycine 50 . As shown in Supplementary  Figs S3 and S4, different members in the GSTU gene family were found to be differentially expressed and the content of glycine were both decreased under the Pb or Cd exposure of radish. Due to the diversity of structural classes of metabolites in plants, there is so far no single technique to identify and quantify all of them 15 . Here, the most critical biological thiols compound involved in GSH metabolism such as GSH, phytochelatins (PCs) and other sulfides could not be successfully detected by the GC-MS. In the previous studies, all the quantification of GSH and PCs were carried out through a targeted-metabolic-profiling analysis by LC-MS method 48,51 . In order to detect a wide range of metabolites, other techniques including LC-MS would be employed to investigate the metabolic profiles in radish upon heavy metal stress, and then more comprehensive profiles would be obtained with the different techniques.
In summary, HM-induced metabolomic responses by GC-MS analysis for radish roots during Pb and Cd stress showed different but overlapping metabolic alterations mainly on sugars, amino acids and organic acids. By integrated analysis of transcript and metabolic profiles, it demonstrate that exposure of radish to Pb or Cd stress resulted in profound biochemical pathway changes, and the integrative gene-to-metabolite networks were revealed in this study. These findings may provide useful information for understanding the molecular mechanisms involved in radish responding to HM-induced stresses, which would further facilitate the effective management of HM contamination in vegetable crops by genetic manipulation.

Methods
Plant materials. Seeds of an inbred line of radish (Raphanus sativus L.), 'NAU-RG' were surface-sterilized for 2 min with 2.5% NaClO followed by washing and soaking in sterile distilled water. After germinating in the dark at 25 °C, the seeds were transferred into culture bowls in a growth chamber (24/16 °C day/night temperature, 14/10 h light/dark). Seedlings with four fully expanded true leaves were transplanted into the culture pots, which were filled with modified half-strength Hoagland nutrient solution, as previously described 24 . After two weeks, the seedlings were treated with either 0, 1,000 mg · L −1 Pb(NO 3 ) 2 , or 400 mg · L −1 CdCl 2 •2.5H 2 O for 72 h. Eight biological replicates from each treatment were sampled and handled according to the reported descriptions 24,25,52 . For each replicate with five seedlings, nearly equal amount of samples from three randomly selected individual plants were pooled. All the harvested radish taproot samples were immediately frozen in liquid nitrogen and then stored at − 80 °C for two weeks before metabolite extraction.
Metabolite extraction from the radish taproots. Metabolites were extracted according to previous studies with some modifications 53,54 . Taproots were immediately plunged into liquid nitrogen to store for further analysis. Briefly, about 50 mg of each frozen samples were transferred into a 2-mL centrifuge tube, and then, 1mL of 100%-methanol (pre-cooled at − 20 °C) was added to each tube. The mixtures were vortexed and ground Scientific RepoRts | 5:18296 | DOI: 10.1038/srep18296 by a 70-Hz grinding mill system (Jinxin Biotech Ltd., Shanghai, China) for 5 min. The homogenate was treated withultra-sonication for 30 min at 70 °C. Subsequently, the tubes were centrifuged at 14,000 g at 4 °C for 10 min. Of the supernatant, 0.4 mL was decanted into a 2-mL screw-top tube, including 10 μ L internal standards (0.02 mg•mL −1 3,4-dichloropheny-lalanine in methanol), 200 μ L chloroform (pre-cooled at − 20 °C) and 400 μ L of demineralized water (Milli Q). Afterwards, the mixtures were vortexed thoroughly and centrifuged for 15 min at 2,200 rcf at 4 °C. Finally, liquids in equal amount of 200 μ L from the aqueous and chloroform layers were transferred into a glass vial for vacuum-dry at room temperature. In addition, Blank samples were also prepared by extraction solution and pooled samples using mixing aliquots from each biological of all the three group taproot samples. Such blank samples were analyzed in each analytic run along with the true samples.
The dried samples were dissolved and derivatized using a two-step procedure involving oximation and silyaltion before injection of GC-MS analysis. First, 30 μ L of methoximation mixture consisting of methoxylamine hydrochloride dissolved in pyridine (20 mg•mL −1 ) were added to the vial, vortexed for 30 s and reacted for 90 min at 37 °C. This was followed by trimethylsilylation with 30 μ L of N,O-bis (trimethylsilyl) trifluoroacetamide (BSTFA) containing 1% trimethylchlorosilane (TMCS), keeping the temperature constant at 70 °C for 60 min. Finally, derived samples were cooled to room temperature before injection.

Gas chromatography-mass spectrometry (GC-MS) analysis. Analysis was performed on an
7890A-5975C GC-MS system (Agilent Technologies, Santa Clara, CA, USA) equipped with an HP-5MS capillary column (30 m × 0.25 mm × 0.25 μ m) (Agilent J &W Scientific, Folsom, CA). All the samples and replicates were continuously injected as one batch in random order to discriminate technical from biological variations. Additionally, the prepared pooled samples were used as quality controls (QCs), which were injected at regular intervals (every ten samples) throughout the analytical run to provide a set of data from which the repeatability can be assessed.
Each 1 μ L aliquot of the derivatized sample solution was injected in splitless into the GC column by a TriPlus autosampler (Thermo Fisher Scientific), at a constant flow rate of helium of 1 mL·min −1 . The temperature of the injector was operated isothermally at 280 °C. The mass spectrometry (MS) transfer line temperature to the quadruple was set to 150 °C and the electron impact (EI) ion source ion temperature was 230 °C. Compound elution settings were 2min at 80 °C isothermal, followed by a 10 °C•min −1 oven temperature gradient to a final 325 °C, and then hold for 6 min at 325 °C. The system is then temperature equilibrated for 1 min at 80 °C before injecting the next sample. Ions were generated by a 70 eV electron beam. The spectra were recorded with a scanning range of 50-550 m • z −1 .
Data analysis and metabolic pathway construction. The '.CDF' formats files obtained from the mass spectra, were imported into the XCMS software which is based on the R software platform (http://cran.r-project. org). XCMS could be employed to be preprocessed in an automatic way including raw signal extracting, data baseline filtering, peak identification, and integration 55,56 . Basically, after alignment with statistic compare component, the '.CSV' file was obtained with three dimension data sets including sample information, retention time and peak intensities. The internal standard was used for data quality control (reproducibility). Internal standards and any known pseudo positive peaks, such as peaks caused by noise, column bleed and BSTFA derivatization procedure, were removed from the data set. Furthermore, the data set was normalized using the sum intensity of the peaks in each sample, which were separately imported into SIMCA-P software package (version 11.0, http://www.umetrics. com/simca). Principal component analysis (PCA) and partial least-squares-discriminant analysis (PLS-DA) were applied to the data after mean-centering and unit variance scaling (UV-scaling). These analyses employed a default seven-fold internal cross validation from which the R 2 and Q 2 (goodness of prediction) values, representing the total explained variance and the model predictability, respectively, were extracted 57,58 .
Metabolites were identified by comparing the mass-to-charge ratios and the abundance of each compound detected against a standard mass chromatogram in the National Institute of Standards and Technology (NIST) database and the Wiley registry of mass spectral library. Peaks with similarity index more than 70% were tentatively identified as metabolites, whereas those having less than 70% similarity index were considered as unknown metabolites. Two-sample t-test statistics was used for the comparison of the metabolite levels to determine their significant differences. Relationships between the significantly altered metabolites were analyzed with R software (KEGGSOAP package). Metabolites among five steps of reactions were considered related ones, and the visualized networks of metabolite-to-metabolite and gene-to-metabolite were created with cytoscape software 59,60 .