Unbiased characterization of genotype-dependent metabolic regulations by metabolomic approach in Arabidopsis thaliana

Metabolic profiles vary between WT, mto1 and tt4 mutants

The metabolome of the aerial regions of individual WT plants and the mto1 and tt4 mutants were analyzed by GC-TOF/MS. The obtained raw data were imported to the MATLAB software for peak alignment and for extraction of spectra using a hierarchical multivariate curve resolution (H-MCR) method as described previously [26, 27]. By this method, we could obtain faster a correct data table with peak areas, including corresponding mass spectra, than by commercial software often used for metabolome analysis [28]. The H-MCR process extracted 518 metabolite peaks and their mass spectra. Using the national institute of standards and technology (NIST02) and our custom software for spectral annotation (Fukushima et al., in preparation), 93 peaks were identified or annotated as known metabolites. In addition, 78 peaks were annotated with mass spectral tags (MSTs) [29, 30]. These peaks were consistently observed though not identified completely.

For a total of 171 annotated peaks, partial least square-discriminate analysis (PLS-DA) was applied as a multivariate statistical analysis. The scatter plot of PLS-DA scores showed clear separation among WT, mto1, and tt4 with a good cross-validation result in accordance with their genetic background (Figure 1). The 2 mutants showed distinct metabolic snapshots that were consistent within each biological replicate of mutant plants.

Subsequently, we searched the metabolites that contributed to this separation of the mutant profiles group from those of WT (discriminative metabolites, hereafter) (Additional File 1). In order to filter the discriminative metabolites further, Welch's t-test was conducted for annotated peaks by interpreting the first weight vector in the PLS-DA, assuming that metabolite concentrations follow the normal distribution in the same genetic background. Between mto1 and WT, 34 peaks differed significantly, whereas 31 peaks were significantly different between tt4 and WT (p < 0.05). The discriminative metabolites for mto1 were mostly amino and organic acids, whereas sugars and the precursors of secondary products were found as discriminative metabolites for tt4. In mto1, remarkable changes in methionine-related compounds such as methionine, homocysteine, and methionine sulfone were observed (see below).

Pronounced metabolite changes are observed in mto1 but fewer changes in tt4

In order to identify the metabolic pathways that were modulated in terms of metabolite levels, significant metabolite changes (p < 0.05) observed in mto1 or tt4 against WT were compared on the metabolic map (Figure 2). Only the peaks annotated as known compounds are projected in Figure 2, while a comprehensive list of all measured peaks and the changes of these levels are shown in Additional File 2.

In mto1, a significant increase was observed in the levels of methionine (34-fold), homocysteine (44-fold), S-adenosyl-methionine (2.3-fold), and methionine sulfone (12-fold) (Figure 2A). Upregulation was also observed in the levels of ornithine, N-acetyl-ornithine, and arginine in the glutamate family. The levels of other amino acids in the aspartate family (aspartate, homoserine, and threonine) were nearly consistent with the reported data [23, 24, 31]. In contrast, there was a significant decrease in the levels of pipecolate (0.4-fold) and uracil (0.5-fold).

Less pronounced changes were observed in tt4 as compared to mto1. The levels of major sugar pools, such as sucrose, glucose, and fructose, showed a significant decrease (Figure 2B). In addition, a few components of the tricarboxylic acid (TCA) cycle and the shikimate pathway were slightly decreased in tt4. The levels of shikimate, sinapate, and tyramine were decreased but that of tyrosine showed no change. The levels of all the amino acids remained unchanged in tt4, suggesting that a lack of flavonoid biosynthesis does not largely affect the metabolism of amino acids.

Metabolite correlations are similar between WT and tt4 but not mto1

In order to obtain insights into the regulation of metabolic network, we focused on the changes in metabolite-to-metabolite correlations in WT and the 2 mutants. Considering the small fluctuations in metabolite levels across individual plants under a given growth condition, metabolite correlations are observed as a representative metabolic profile of a given genotype [12, 14, 32, 33]. The correlation matrices of the WT plant and the 2 mutants were visualized by heat maps (Figure 3). The numbers of significantly correlated metabolite pairs in each genotype plant (1250 in WT, 305 in mto1, and 3751 in tt4; p < 0.001) in Figure 3 indicated that mto1 lost metabolite correlations and tt4 strengthened correlations when compared with the WT plant, though WT and 2 mutants showed similar growth and phenotypes in the experimental condition, which was strictly controlled as much as possible. In WT and tt4, fatty acids had a trend of negative correlation with all the other components. In contrast, fatty acids and steroids showed positive correlations with each other. For the sugar groups in WT and tt4, fructose and glucose had a tendency of negative correlations with the other components; however, fructose-glucose correlated in a highly positive manner. Although mto1 had lesser metabolite correlations as compared to the other 2 genotypes, characteristic correlations were found only in mto1, e.g., linolenate (or linoleate)-steroids, linolenate (or linoleate)-nicotianamine, and glutamate-glutamine. The correlations specific to tt4 were also observed, such as malate-aromatic amino acids (see below).

The correlations obtained in each genotype were filtered by the statistical significance of the correlation coefficient, and these significant pairs were compared among the 3 genotypes (Figure 4). The correlation profile of mto1 differed largely from those of the remaining 2 genotypes. Approximately 8.3% of the significant correlation pairs were common among WT, mto1, and tt4; 88.4% (80.1% plus 8.3%) of the pairs in WT were conserved in tt4; only 9.1% (0.8% plus 8.3%) of the pairs were conserved between WT and mto1. These results suggested that the deregulated overproduction of methionine in mto1 strongly affected the stability of the cellular metabolite network, whereas the lack of flavonoids in tt4 resulted in the strengthening of metabolite correlation networks.

Malate and aromatic compounds in shikimate pathways are tightly coregulated in tt4

Among the metabolite correlation pairs significant in tt4, malate and sinapate exhibited a much tighter correlation (correlation coefficient: r _Met = 0.89) than did WT (r _Met = 0.50) and mto1 (r _Met = 0.06) although these 2 metabolites are located at a distance in the metabolic map. Table 1 shows the correlation of malate or sinapate with the other metabolites belonging to the TCA cycle and shikimate pathways. Malate- citrate and malate- fumarate, which are components in the TCA cycle, exhibited significant correlations both in WT and tt4. On the other hand, malate and the aromatic compounds in the shikimate pathway (phenylalanine, tyramine, tyrosine, and sinapate) were correlated only in tt4 but not in WT (Table 1 and Figure 5) although the levels of these metabolites in tt4 were maintained at levels similar to those observed in WT (see Additional File 2). Among the compounds concerned with these unique correlations enforced in tt4, malate and sinapate are believed to be direct precursors for the formation of sinapoylmalate acting as a UV-B protectant [34, 35]. Interestingly, the level of sinapoylmalate itself did not change in tt4 as compared with that in WT (Additional File 3). These results suggested that the correlations between malate- aromatic compounds were intensified in tt4 without any significant changes in the pools of metabolites concerned (Figure 5). The tight correlation between malate and sinapate may imply a connection with the role of sinapoylmalate against UV-B stress in the tt4 mutant, in which the UV-protective flavonoids are missing (see Discussion).

Metabolite X	Metabolite Y	r Met in WT	r Met in mto1	r Met in tt4	Pathway
Malate	Maleate	0.80	0.46	0.91	TCA
	Succinate	0.70	-0.06	0.89	TCA
	Fumarate	0.83	0.07	0.90	TCA
	Citrate	0.92	0.79	0.87	TCA
	Phenylalanine	0.45	-0.01	0.89	Shikimate
	Shikimate	0.81	0.57	0.92	Shikimate
	Tyramine	0.44	0.10	0.94	Shikimate
	Tyrosine	0.48	0.48	0.92	Shikimate
	Sinapate	0.50	0.06	0.89	Shikimate
Sinapate	Maleate	0.74	0.59	0.90	TCA
	Succinate	0.85	0.74	0.86	TCA
	Fumarate	0.75	0.87	0.90	TCA
	Citrate	0.36	0.09	0.84	TCA
	Malate	0.50	0.06	0.89	Shikimate
	Phenylalanine	0.74	0.33	0.86	Shikimate
	Shikimate	0.72	0.38	0.88	Shikimate
	Tyramine	0.80	0.11	0.91	Shikimate
	Tyrosine	0.71	0.29	0.89	Shikimate

Boldface characters represent the unique correlations observed only in tt4. The unique correlations in tt4 are defined as r _Met > 0.85 in tt4 and as r _Met < 0.80 in WT and mto1. Abbreviations: TCA, tricarboxylic acid.

Few metabolite correlations are conserved among WT, mto1, and tt4 mutants

Figure 6A presents the metabolite correlation network in WT. The metabolite correlation pairs conserved among WT, mto1, and tt4 are listed in Table 2 (only for identified metabolites) and are shown in Figure 6B (all the detected peaks). These correlation pairs conserved in the 3 genotypes are categorized into the following 3 types with respect to biosynthetic relations: (1) metabolite pairs placed close to each other in a metabolic pathway for the concerned reactions, e.g., the pairs of glucose-6-phosphate- fructose-6-phosphate, fructose- glucose, and the component pairs of the TCA cycle; (2) metabolite pairs possessing a similar chemical skeleton derived from the same biosynthetic pathway, e.g., the pairs of cholesterol, campesterol, and beta-sitosterol; (3) metabolite pairs located at a distance in a metabolic map, e.g., the pairs of valine- beta-alanine and serine- glutamate. Both valine and beta-alanine belong to the pantothenate (vitamin B₅) synthesis process, although current knowledge regarding their regulatory relationship in this pathway in plants is limited. Serine and glutamate are involved in many metabolic reactions in various biosynthetic pathways and are directly connected by serine aminotransferase. These conserved metabolite correlations regardless of the mutations in mto1 and tt4 are presumed to reflect the fundamental metabolic regulations underlying plant metabolic systems. It is also notable that the network clusters in Figure 6 adequately represent the metabolic pathways, e.g., sucrose degradation in (1), steroids in (2), and the TCA cycle and amino acids in (3), indicating that the correlation network forms a modular structure, at least partly, according to the metabolic pathways.

A new score for metabolites can select metabolomic markers of a mutant

In order to represent the metabolic characters of plants with different genotypes particularly in terms of metabolite correlation besides simple comparison of changes in metabolite levels, a new index for a metabolite representing the change of metabolite correlation was necessary. We therefore defined a new scoring method for metabolites that represents a diverse correlation in a genotype. This new criterion "z- score of metabolite correlation (ZMC)" represents the extent of uniqueness over given genotypes in terms of metabolite correlation, i.e., if metabolite correlation coefficient of metabolite X in a given genotype decreases markedly as compared to those in other genotypes, the ZMC of X in this genotype increases.

It should also be noted that the simple comparison of these metabolite levels listed by ZMC did not show significant differences except in the case of methionine and methionine sulfone. Thus, combinatory analyses of metabolite levels and metabolite correlation, such as ZMC, are important to obtain a clear insight into the manner in which the system is regulated and to generate a new hypothesis for characterizing metabolic systems in plants with different genotypes.

Projection of metabolite correlation to gene-expression correlation gives clues on how metabolic networks are regulated

Since the conserved metabolite correlations across the different genotypes are presumed to characterize fundamental networks regardless of mutations, the investigation of gene expression in the metabolism of these metabolite pairs may lead to deeper insight into network regulation. The nature of metabolic networks represented as metabolite correlation might be explained, at least partly, by the correlation of gene expressions, and new findings regarding the gene-expression network could be obtained from the correlations of conserved metabolite.

Figure 7 illustrates the flow chart of the mining of gene-expression correlation based on metabolite correlation pairs. Given metabolites X and Y as a correlating pair, we initially searched for genes coding for enzymes directly involved in the metabolic reactions of metabolite X or Y acting as a reactant or as a product in AraCyc metabolic map [36]. Subsequently, coexpression analysis among extracted enzyme-coding genes was performed by using the ATTED-II coexpression database of Arabidopsis transcriptome [22]. Coexpressed patterns (r _Exp ≥ 0.6) were observed in 49 gene-expression pairs across the whole transcriptome database (listed in Additional File 4). In order to exclude the pairs whose apparent coexpression pairs are simply due to their narrow tissue-specific expression patterns, we searched for Gene Atlas in Genevestigator [21] for the expression of condition specificity. Of these 49 enzyme-coding gene pairs, 6 coexpression pairs were indeed simply explained by tissue-specific expression patterns (for example, pollen and stamen), while the remaining 43 pairs of coexpression were most likely ascribed to their intrinsic coexpression natures, presumably causing the metabolite correlations. These gene pairs would prove interesting for further investigation to elucidate the bases of metabolic networks such as finding the transcription factors coexpressing with these enzyme-coding genes (as also shown in Figure 7). However, if we calculate the proportion of coexpressed enzyme-coding gene pairs over all the possible combinations of genes involved, only a very small fraction of combinations (0.2–3.9%) showed coexpression (Additional File 5) indicating the presence of multiple layers of regulation of the metabolic network not restricted to the level of gene expression. For further results see Additional file 6.

Single deregulation of methionine synthesis disturbs the entire primary metabolism

For mto1, a significant increase in the level of homocysteine, methionine, and S-adenosyl-methionine was observed as expected from the overaccumulation of methionine; unexpectedly, there was a significant decrease in the levels of pipecolate and uracil (Figure 2A). Pipecolate is known to be a component of the lysine degradation pathway in a variety of organisms including mammals, fungi, and plants [37–39]. Uracil is derived from the aspartic acid family connected to methionine. In addition, there was a significant increase in the levels of arginine and ornithine. Arginine and ornithine belong to the glutamate family and are precursors of putrescine, which is an important metabolite in polyamine biosynthesis [40]. The results of ZMC analysis also indicated spermidine as an index metabolite of mto1. These data suggest that the overaccumulation of methionine affects not only methionine metabolism but also the regulation of the glutamate family leading to the regulation of polyamines. With regard to the correlations of metabolites belonging to the aspartate, glutamate, and polyamine biosynthetic pathways, most of the metabolites were highly conserved in both WT and tt4 but not in mto1. These results indicate that the perturbation in methionine metabolism causes a strong impact and a reconfiguration in the cellular nitrogen status by finely adjusting the metabolic networks of the nitrogenous metabolites.

Backup system is reconfigured for protection against UV stress in the tt4 mutant

The metabolite pairs of the malate- shikimate pathways, in particular, the malate- sinapate pair, exhibited enhanced correlations in tt4 although the levels of these compounds remain unchanged (Figure 5). Sinapoylmalate, a major metabolite of sinapoylesters, is believed to play a protective role against UV-B in Arabidopsis [41]. The correlation observed between malate and sinapate in tt4, which is not observed in WT, implies the reconfiguration of metabolic networks lacking flavonoids in tt4 to efficiently produce an alternative defense compound, sinapoylmalate, thereby replacing the UV-protectant flavonoids. It is worthwhile to note that this correlation was observed even under normal conditions of light, in which no elevated accumulations of sinapoylmalate and its precursors were observed; this suggests the presence of a yet unidentified mechanism of cellular preparation for an adaptive response to UV stress by the reconfiguration of the metabolite networks in Arabidopsis. When we performed the UV-B stress experiment on WT and tt4, the levels of sinapoylmalate and sinapoylglucose were increased more in tt4 than in WT under UV irradiation (unpublished results). More detailed statistical and correlation analyses for sinapoylmalate, sinapoylglucose, and other aromatic compounds are in progress.

A novel strategy for the prediction of genes responsible for metabolic networks by combined correlation analysis of metabolite accumulation and gene expression

There are a few reports available in literature concerning the correlations of metabolite accumulation, enzyme activities, and gene expression. The metabolite correlation, mainly in primary metabolites, is presumed to represent the direct and indirect interconnections in the metabolic pathway, even in the absence of causality between the correlation profiles and the pathway [18]. Non-synchronized patterns of transcripts and enzymes have been observed in WT and the starch-free mutant of Arabidopsis [42, 43]. The significant metabolite-to-gene correlation is likely to reflect the effects of metabolic regulation that occur from metabolites to gene expression, rather than those from gene expression to metabolites [13]. The origin of high metabolite correlation is suggested to be derived either from the strong control of a single common enzyme between metabolites or from the increase in the level of a single enzyme as compared to the others [44].

Several studies indicate the power of integrated analysis of the metabolome and transcriptome [9–11], in particular, the feasibility of coexpression analysis by public transcriptome datasets [45, 46], for gene characterization. Thus, we conceived that the basis of metabolic networks conserved across plants of different genotypes can be explained, at least partly, by the correlation of gene expression, and a better understanding of gene function can be obtained in the context of metabolic networks.

When we searched for the coexpression of enzyme-coding genes directly connected to metabolite pairs whose correlation is conserved, only a small fraction of genes were coexpressed. Nevertheless, we could filter the coexpressed enzyme-coding gene pairs (Additional File 4) that are presumably important for the regulation of metabolic networks. With regard to the serine-glutamate pair, serine hydroxymethyltransferase 1 (SHM1) (At4g37930) showed a strong correlation with glutamate:glyoxylate aminotransferase 1 (GGT1) (At1g23310). Both enzymes are involved in the photorespiration/metabolic salvage pathway, suggesting a coordinate regulation of this pathway at both the gene expression and metabolite levels. In addition, a number of genes of unknown function are also listed. It is intriguing to deduce the function of these genes by accounting for their correlation with metabolite pairs by further study.

Furthermore, we extended the coexpression analysis to the transcription-factor genes with the enzyme-coding genes (Figure 7). This search extracted CONSTANTS-LIKE3 (COL3) (At2g24790) and 3 unknown transcription factor genes (At1g68520, At4g36540, and At4g14540) from the metabolite pair of serine-glutamate. COL3 was previously identified as a constitutive photomorphogenetic 1 (COP1)-interacting protein and was suggested to be a positive regulator of light signaling [47]. Three unknown genes exhibited many connections with several enzyme-coding genes (SHM1, GGT1, glutamine synthetase 2 [GS2], and glutamate-1-semialdehyde 2,1-aminomutase 2 [GSA2]) involved in photorespiration/metabolic salvage pathway and nitrate assimilation. These unknown transcription-factor genes are of interest for further characterization because these may play an important role in the metabolic network of serine-glutamate. In addition to the serine-glutamate pairs, the other metabolite correlation pairs led to the suggestion of several transcription-factor genes through enzyme-coding gene coexpression in the same manner. This strategy can be a novel way of delimiting genes involved in metabolic networks.

Reagents

All chemicals excluding isotope reference compounds and reagents for silylation were purchased from Sigma Aldrich (Tokyo, Japan). The 6 stable isotope compounds ([¹³C₅]-proline, [²H₄]-succinic acid, [²H₆]-2-hydroxybenzoic acid, [¹³C₃]-myristic acid, [¹³C₁₂]-sucrose, [²H₇]-cholesterol) were purchased from Cambridge Isotope Laboratories (Andover, MA, USA), and [¹³C₅, ¹⁵N] glutamic acid and [¹³C₆]-glucose were obtained from Spectra Stable Isotopes (Columbia, Maryland, USA). [²H₄]-1,4-diaminobutane was from C/D/N ISOTOPES (Pointe-Claire, Quebec, Canada). [¹³C₄]-hexadecanoic acid was from Icon (Mt. Marion, NY, USA). The reagent for trimethylsilylation, N-methyl-N-trimethylsilyltrifluoroacetamide (MSTFA) with 1% of trimethylchlorosilane (TMCS), was purchased from Pierce (Rockford, IL, USA).

Plant material and harvest

Wild-type A.thaliana plants accession Columbia (Col-0) and the mutants mto1 [23] and tt4 [48] of Col-0 background were obtained from Dr. Naito, Hokkaido University, and Dr. Kitamura, Japan Atomic Energy Research Institute, respectively. The sterilized seeds were stratified at 5°C for 2 d, and were successively sown on Murashige and Skoog (MS) medium containing 1% sucrose. Plants were cultivated in controlled growth chambers at 22°C in 16-h light and 8-h dark conditions for 18 d. The aerial regions were harvested with 20 different biological replicates, 6 h after the onset of the bright phase. Among these harvested plants, the plants that had fresh weight over 13 mg were used for GC-TOF/MS analysis (WT, n = 17; mto1, n = 13; and tt4, n = 20). For liquid chromatography-quadrupole-time-of-flight/mass spectrometry (LC-Q-TOF/MS) analysis, 3 biological replicates were prepared. All the plant materials were frozen immediately in liquid nitrogen to quench the enzymatic activity.

Extraction and sample preparation for GC-TOF/MS analysis

Each sample was extracted with a concentration of 25 mg fresh weight of tissues per μl extraction medium (methanol/chloroform/water [3:1:1 v/v/v]) containing 10 stable isotope reference compounds using a Retsch mixer mill MM310 at a frequency of 30 Hz for 3 min at 4°C. Each isotope was adjusted to a final concentration of 22.5 ng per μl injection [49, 50]. After centrifugation for 5 min at 15,100 × g, 400 μl of the supernatant was used for further analysis. The extracts were evaporated to dryness in a Savant SPD2010 SpeedVac Concentrator (Thermo Electron Corporation, Waltham, MA, USA).

For methoximation, 20 μl of methoxyamine hydrochloride (20 mg/ml in pyridine) was added to the sample. After 30 h of derivatization at room temperature, the sample was trimethylsilylated for 1 h using 20 μl of MSTFA with 1% TMCS at 37°C with shaking. Twenty μl of n-heptane was added following silylation. All the derivatization steps were performed in the vacuum glove box VSC-100 (Sanplatec, Japan) filled with 99.9995% (G3 grade) of dry nitrogen. The analysis of metabolites by GC-TOF/MS was performed as described previously [28].

LC-Q-TOF/MS analysis

Frozen leaves were homogenized in 5-μl extraction solvent (methanol/H₂O [4:1 v/v]) per mg fresh weight of tissues by using a mixer mill at a frequency of 20 Hz for 3 min at 4°C. After centrifugation at 12,000 × g, the cell debris was discarded, and the extracts were centrifuged again. These supernatants were immediately used for flavonoid analysis. For flavonoid profiling, Waters Acquity UPLC™ system (Waters Co., Massachusetts, USA) fitted with a Q-ToF Premier mass spectrometer (Micromass MS Technologies, Manchester, UK) was used. Ultra-performance liquid chromatography (UPLC) was carried out on a UPLC™ BEH C₁₈ column (100-mm length × 2.1-mm inner diameter, 1.7-μm particles, Waters Co.) at a flow rate of 0.5 ml/min at 35°C. The elution gradient comprised solvent A (0.1% trifluoroacetic acid in H₂O) and solvent B (0.1% trifluoroacetic acid in acetonitrile) and the elution profile – 0 min, 100% A; 5 min, 11% B; 20 min, 13% B; 24 min, 100% B – using linear gradients in between the time points. A photodiode array (PDA) detector was used for the detection of UV-visible absorption in the range of 210–500 nm. The TOF mass analyzer was used for the detection of flavonoid glycosides [M+H]⁺ and fragment ions peak in a positive ion mode scan. The desolvation temperature was 450°C with a nitrogen gas flow rate of 600 l/h, capillary spray at 3.2 kV, source temperature at 150°C, and cone voltage at 35 V.

The peaks in the plant extracts were identified based on retention times, UV visible absorption spectra, and mass fragmentation by tandem MS analysis as reported [45]. The amounts of each kaempferol glycoside and sinapoyl derivative were calculated using kaempferol (at 340 nm) or sinapic acid (at 340 nm) as a standard, respectively.

Processing of GC-TOF/MS data

Nonprocessed MS data from GC-TOF/MS analysis were exported in NetCDF format to MATLAB 6.5 (Mathworks, Natick, MA, USA), where all data-pretreatment procedures, such as smoothing, alignment, time-window setting, and H-MCR, were carried out [27]. The resolved MS spectra were matched against reference mass spectra using the NIST mass spectral search program for the NIST/EPA/NIH mass spectral library (version 2.0) and our custom software for peak annotation written in JAVA. Peaks were identified or annotated based on retention indices (RIs) and the reference mass spectra comparison to the Golm Metabolome Database (GMD) [30] and our in-house spectral library. The metabolites were identified by comparison with RIs from the library databases (GMD and our own library) and with those of authentic standards, and the metabolites were defined as annotated metabolites on comparison with mass spectra and RIs from these two libraries. The amount of S-adenosyl-methionine was calculated by the sum of the mass numbers at m/z 188 and 236 using Leco ChromaTOF software version 2.32 (LECO, St. Joseph, MI, USA) since this compound was not adequately detected by H-MCR.

Statistical data analysis for metabolic profiling and metabolite correlation analysis

After peak annotation, the obtained data matrix (observations: biological replicates, variables: 171 annotated peaks) was used for statistical analyses including multivariate analysis and correlation analysis as log₁₀-transformed data. Multivariate analysis was performed with SIMCA-P 11.0 software (Umetrics AB, Umeå, Sweden). Two types of multivariate analyses were applied: principal component analysis (PCA) as an unsupervised analysis and partial least square-discriminate analysis (PLS-DA) as a supervised analysis (WT, n = 17; mto1, n = 16 [13 samples from biological replicates and 3 samples from analytical replicates]; and tt4, n = 20). The PCA was applied to check the sample distribution before performing PLS-DA. In order to narrow down the loadings that make a contribution to separate the groups in PLS-DA, the loading of principal component 1 in each model was chosen according to the value of the first weight vector (w* 1) together with the 99% confidence intervals calculated using jack-knifing [51, 52]. Statistical tests were performed using R package, which is an open source statistical computing environment [54]. Simple comparisons of the means of the obtained peak areas were performed by Welch's t- test. A difference with p < 0.05 was considered to be significant. Pair-wise metabolite- metabolite correlations were calculated by Pearson's correlation coefficient (r _Met ) test using 171 peaks. The level of significance was determined as follows: |r _Met | ≥ 0.8 (p < 1.2 × 10^-4, n = 17) for WT, |r _Met | ≥ 0.85 (p < 2.3 × 10^-4, n = 13) for mto1, and |r _Met | ≥ 0.75 (p < 1.4 × 10^-4, n = 20) for tt4. The network of pair-wise metabolite correlation is drawn using Pajek [55].

ZMC was calculated as follows: (1) For a metabolite in a given genotype (e.g., mto1), all the significant metabolite correlation pairs (e.g. r _Met ≥ 0.8) which contained the metabolite as a pair in 2 genotypes (e.g., WT and tt4) were chosen; (2) For the pairs selected in step (1), the average correlation coefficient in the 2 genotypes (e.g., WT and tt4) was calculated; (3) The difference in the correlation coefficient between a given genotype and the other 2 genotypes was obtained by subtracting the average coefficient in 2 genotypes from the coefficient of a given genotype; (4) The z-score of each relative correlation obtained was calculated and defined as a ZMC. Finally, the ZMCs across all pairs were computed. Given that ZMC follows normal distribution, the ZMC of a metabolite is regarded as an index representing the characteristic metabolotype in terms of metabolite correlation in a given genotype.

Coexpression analysis of genes using publicly available data source

In order to understand the relationship between metabolite-metabolite correlations and the coexpression patterns of the genes encoding collateral enzymes for each concerned metabolite, publicly available data that were obtained from AraCyc database version 3.5 [36] as a text format dump file "AraCyc Pathways" that listed accessions for 262 different pathways associated with 1885 genes (aracyc_dump_20070213) and another file "AraCyc Compounds" that contained all compounds found in the AraCyc pathways (aracyc_compounds_20070213) were used. When metabolite X (or metabolite Y) reacts with an enzyme as a substrate or a product, the enzyme was defined as a "reaction" enzyme. The coexpression analysis was performed between the enzyme-coding genes for metabolite X and metabolite Y. In order to search for coexpressed genes, we used a tool named "Correlated Gene Search" from the website [56] released by our institute (to be published elsewhere). The expression data were produced using Affymetrix GeneChip by AtGenExpress at the Arabidopsis information resource (TAIR), and the correlation data has been released in the ATTED-II database [22]. The search conditions were as follows: the type of matrix was "all dataset ver. 3" since the source of GeneChip data were the 58 experiments including 1,388 arrays; threshold = 0.6, method = "interconnection of sets", display limit = 99999 (maximum). The distribution of the correlation coefficient for coexpressed genes via ATTED-II followed the normal distribution (the standard deviation was equal to 0.18). The result suggested that a threshold above 0.6 could correspond with 3 standard deviations.

For a coexpression analysis of the transcription factor and coexpressed enzyme-coding gene pairs across the reactions for each metabolite, a list of 1968 transcription-factor genes was made from RIKEN Arabidopsis Transcription Factor database (RARTF) [53]. Using AGI code numbers in the list, coexpression analysis was performed using the tool "Correlated Gene Search" from the website [56]. The organ-specificity of the expression of the coexpressed genes was investigated in silico by using Genevestigator [21].