Metabolic network fluxes in heterotrophic Arabidopsis cells: stability of the flux distribution under different oxygenation conditions.

Steady-state labeling experiments with [1-(13)C]Glc were used to measure multiple metabolic fluxes through the pathways of central metabolism in a heterotrophic cell suspension culture of Arabidopsis (Arabidopsis thaliana). The protocol was based on in silico modeling to establish the optimal labeled precursor, validation of the isotopic and metabolic steady state, extensive nuclear magnetic resonance analysis of the redistribution of label into soluble metabolites, starch, and protein, and a comprehensive set of biomass measurements. Following a simple modification of the cell culture procedure, cells were grown at two oxygen concentrations, and flux maps of central metabolism were constructed on the basis of replicated experiments and rigorous statistical analysis. Increased growth rate at the higher O(2) concentration was associated with an increase in fluxes throughout the network, and this was achieved without any significant change in relative fluxes despite differences in the metabolite profile of organic acids, amino acids, and carbohydrates. The balance between biosynthesis and respiration within the tricarboxylic acid cycle was unchanged, with 38% +/- 5% of carbon entering used for biosynthesis under standard O(2) conditions and 33% +/- 2% under elevated O(2). These results add to the emerging picture of the stability of the central metabolic network and its capacity to respond to physiological perturbations with the minimum of rearrangement. The lack of correlation between the change in metabolite profile, which implied significant disruption of the metabolic network following the alteration in the oxygen supply, and the unchanging flux distribution highlights a potential difficulty in the interpretation of metabolomic data.


INTRODUCTION
Although the complexity and plasticity of the metabolic network in plants allows them to adapt to fluctuating environmental conditions, the same properties also present a significant obstacle to metabolic engineering (Carrari et al., 2003a;Kruger and Ratcliffe, 2008;Sweetlove et al., 2008).
The problem is particularly acute in primary metabolism, where there have been numerous instances of unsuccessful engineering, and reflects the current incomplete understanding of the way in which metabolic networks respond to environmental and genetic perturbations. Fluxes of central carbon metabolism are part of the missing information (Sweetlove et al., 2003), and although they are necessarily related to enzyme abundances, metabolite concentrations and transcriptional responses (Carrari et al., 2006;Junker et al., 2007) their reliable prediction from the available data remains a non-trivial task (Sweetlove and Fernie, 2005). For this reason the development and application of techniques for the measurement of flux in plants has become an important area of research (Schwender et al., 2004a;Fernie et al., 2005;Ratcliffe and Shachar-Hill, 2006). Steady state metabolic flux analysis (MFA) has the capacity to resolve parallel, cyclic and reversible fluxes making it a useful technique for quantifying metabolic fluxes, and investigating the factors that control them, in plants (Roscher et al., 2000;Ratcliffe and Shachar-Hill, 2006). MFA studies have revealed novel aspects of plant metabolism, as well as providing the first measurements of many fluxes in vivo and independently verifying previous research on plant metabolism (Schwender et al., 2004a). For example, recent work on Brassica napus embryos highlighted the importance of Rubisco in refixation of CO 2 for improved conversion of photosynthate to seed storage compounds (Schwender et al., 2004b); while work on sunflower embryos (Alonso et al., 2007a) provided evidence that pyruvate uptake is not the dominant route for the provision of precursors for plastidic fatty acid synthesis, in agreement with recent work on Arabidopsis (Andre et al., 2007). Recent refinements to MFA theory (Sriram and Shanks, 2004a;Ghosh et al., 2006;Kruger et al., 2007a;Libourel et al., 2007) and the development of specialised software (Wiechert et al., 2001;Sriram et al., 2004) have been accompanied by a rapid increase in the number of species investigated using MFA (Schwender, 2008). However, the approach has yet to be applied to Arabidopsis thaliana and this paper describes an MFA study of respiratory and biosynthetic carbon metabolism in a heterotrophic Arabidopsis cell suspension.
A key question in the regulation of central carbon metabolism is how the simultaneous demands of catabolic respiratory metabolism and anabolic biosynthetic metabolism are managed. The tricarboxylic acid (TCA) cycle is central to both processes, generating reducing equivalents for the mitochondrial electron transfer chain, and providing precursors for several biosynthetic pathways (Fernie et al., 2004). Recent studies have helped to demonstrate the degree to which the TCA cycle is connected to other metabolic processes: for example, antisense knockdown of malate dehydrogenase (Nunes-Nesi et al., 2005), and fumarase (Nunes-Nesi et al., 2007) and a mutation within the aconitase gene (Carrari et al., 2003b) all had marked effects on photosynthetic performance in tomato, while the effects of oxidative stress on the TCA cycle propagated throughout the metabolic network in A. thaliana (Baxter et al., 2007). MFA has demonstrated that the TCA operates in different flux modes in different systems. For example, in Brassica napus embryos photosynthetic ATP production alleviates the need for any cyclic TCA cycle flux (Schwender et al., 2006), while in other systems there appears to be more conventional operation of the TCA cycle (Rontein et al., 2002, Alonso et al., 2007b with both significant respiratory and biosynthetic flux. MFA can also be used to investigate the response of the TCA cycle and its associated pathways to different growth conditions. Recent work revealed that Brassica napus embryos respond to a switch from an organic to an inorganic nitrogen source by increasing anaplerotic flux through phosphoenolpyruvate carboxylase (PEPC), thereby replacing the carbon removed from the TCA cycle for ammonium assimilation (Junker et al., 2007). In soybean cotyledons relative flux through PEPC was reduced by increased growth temperature (Iyer et al., 2008).
Thus, a quantitative picture is emerging of TCA cycle fluxes, and the extent to which they vary depending on the need to generate precursors for biosynthesis and reductant for ATP synthesis.
However, the factors that may control the rates of biosynthesis and respiration, and the balance between these two competing processes, have not been tested systematically. Accordingly we varied the concentration of O 2 in the medium of an A. thaliana cell suspension culture with the aim of perturbing the operation of the TCA cycle. The effect of this manipulation on the flux map of central metabolism, and in particular the effect on the balance between respiratory and biosynthetic fluxes, was quantified using steady state MFA.

Construction and refinement of a metabolic model
The successful application of MFA requires the construction of a model that not only accurately reflects metabolism within the experimental system, but which can also be solved with the quantity and quality of data that are likely to be obtainable. To this end we constructed an initial model of central carbon metabolism in heterotrophic A. thaliana cell suspension cultures using a format compatible with the steady state MFA software 13C-FLUX (Wiechert et al., 2001). The structure of the network was based on information from the literature, principally other MFA studies of heterotrophic metabolism in plants (Schwender et al., 2006;Alonso et al., 2007a;Sriram et al., 2007), and from metabolic databases (Schomburg et al., 2004;Zhang et al., 2005). Reactions primarily associated with photosynthetic metabolism (Calvin cycle, photorespiration) and seed germination (glyoxylate cycle) were not included since there is no evidence that they occur to any significant extent in dark grown Arabidopsis cell suspensions. Carbon transitions within the network were derived from the primary literature and standard biochemical textbooks. Accurate representation of carbon transitions around the TCA cycle in the model was confirmed by comparing data from a 13 C NMR analysis of the labelling of organic and amino acids in methanolic extracts of heterotrophic cell suspensions fed with [2,3-13 C]succinate for 6 or 18 h with the pattern of label distribution predicted by the model (data not shown).
The complexity of the network that can be analysed is partly determined by the extent to which the redistribution of the 13 C-label can be quantified after labelling to isotopic steady state. We used the statistical analysis component of 13C-FLUX, EstimateStat, to predict errors on optimised flux estimates for different network configurations and 13 C-labelled precursors, and hence refined the initial model to the point where it could be solved with the data obtainable from a steady state labelling experiment. By this method, structurally non-identifiable fluxes, i.e. those that can take any value without impacting on the observed label distribution (Wiechert et al., 2001), were identified and removed from the network. For example, flux from 2-oxoglutarate to succinate, whether occurring via α -ketoglutarate dehydrogenase or via the GABA shunt, was considered as a single flux since these parallel pathways produce identical distributions of label in their product, and sucrose cycling was excluded from the model for the reasons discussed elsewhere (Kruger et al., 2007a). Simplifications were also introduced where the accuracy of flux estimates was predicted to be poor, providing that such simplifications would not prevent conclusions being drawn about the function of the TCA cycle. The refined network is illustrated in Fig. S1.

Determination of optimal 13 C-labelled precursor
EstimateStat was also used to predict the optimal precursor for estimating TCA cycle fluxes in the refined model. Figure 1 suggests that [1-13 C]glucose provides the best estimates of flux through the TCA cycle and elsewhere in the network. [1-13 C]glucose and [U-13 C 6 ]glucose performed similarly well for TCA cycle fluxes, with predicted errors less than the magnitude of the flux estimates, but this analysis suggested that [1-13 C]glucose provided the more accurate estimates of flux elsewhere in the network. [2-13 C]glucose appeared to offer no advantage over [1-13 C] or [U-13 C 6 ]glucose and is significantly more expensive. The analysis was repeated using several different models with more or less explicitly defined subcellular compartmentation of glycolysis and the oxidative pentose phosphate pathway (data not shown). Though the degree of compartmentation greatly affected the predicted relative errors, the qualitative finding that [1-13 C]glucose would provide the smallest relative errors for the majority of fluxes, remained the same. Post hoc analysis of the final model at the end of the investigation (data not shown), confirmed that labelling with 100% [1-13 C]glucose provided the most reliable estimates of flux.

Elevated O 2 conditions perturb cell suspension culture metabolism
As oxygen is required to support respiration, the availability of oxygen to the cell suspension cultures might be expected to influence overall rates of metabolism and/or the partitioning of carbon entering the TCA cycle between respiration and biosynthesis. To test this hypothesis we established a system for culturing cells at elevated O 2 concentration by replacing the aluminium foil used to seal the flasks with Miracloth. After 5 d of growth there was consistently more oxygen dissolved in the medium of cultures covered with Miracloth (elevated O 2 ) than those covered with foil (standard O 2 ): in a representative experiment the oxygen concentration of the culture medium was 161.5 ± 12.5 µM for elevated O 2 cells and 76.0 ± 2.5 µM for standard O 2 cells, both of which are lower than the 270 µM expected for air-saturated water at 21 o C (Truesdale and Downing, 1954). This difference was sufficient to cause significant increases in the abundance of amino acids (glutamate, GABA and alanine) and sugars (glucose) and decreases in the abundance of organic acids (malate, succinate, citrate and fumarate) under elevated O 2 conditions (Fig. 2). In addition, the rates of biomass accumulation and glucose consumption were both greater under elevated O 2 conditions; glucose consumption increased from 214 ± 45 mg d -1 flask -1 under standard conditions to 335 ± 26 mg d -1 flask -1 under elevated conditions, whilst biomass accumulation increased from 174 ± 17 mg d -1 flask -1 to 250 ± 19 mg d -1 flask -1 . However, the relative biomass composition of the cell suspensions after 5 d growth (percentage of biomass that consisted of starch, cell wall, protein, lipids) was the same under both conditions (Table 1). No significant differences could be detected in the composition of the growth medium using 1 H NMR (data not shown). These results suggest that a rearrangement of the metabolic network leads to increased biosynthetic fluxes and a switch from accumulation of organic acids to accumulation of amino acids in response to elevated O 2 . Such a rearrangement could be in the vicinity of the TCA cycle, given its involvement in both respiratory and biosynthetic processes, and this hypothesis was tested using MFA.

Validation of isotopic and metabolic steady state
Steady state flux analysis requires isotopomer abundances to be measured when the system is at isotopic and metabolic steady state, i.e. when metabolic fluxes are constant and when the distribution of label throughout the network has stabilised (Ratcliffe and Shachar-Hill, 2006). To confirm that isotopic steady state was reached after 5 d growth under both elevated and standard O 2 conditions, cells were fed with [U-13 C 6 ]glucose (6% of total glucose supplied) for 4.5 or 5 d and then the distribution of 13 C within soluble metabolites was analysed using 1D 13 C NMR. Cumomer abundances as a percentage of total labelling of a metabolite were calculated and compared between 4.5 and 5 d. There was no significant variation in abundance for the majority of measurements, and although several measurements showed significant variation these reflected changes in abundance of less than 2.5%, strongly suggesting that these metabolites had reached isotopic steady state after 5 d (Fig. 3A). Comparison of percentage cumomer abundances between soluble and protein derived amino acids after labelling with [1-13 C]glucose also revealed few significant differences (Fig. 3B). Such differences as were present reflected changes in abundance of less than 3.9%, suggesting that turnover of protein is sufficiently rapid for it to be labelled to isotopic steady state (Roscher et al, 2000). Thus it was concluded that the cells were at isotopic steady state after 5 d growth.
The existence of an isotopic steady state suggests that any changes in metabolic fluxes over time must be relatively slow (Roscher et al, 2000). Stronger evidence for a metabolic steady state in our system is provided by the fact that the abundance of most soluble metabolites did not change from 4.5 d to 5.5 d (Fig. 2). In particular there were no significant changes in abundance between 5 and 5.5 d growth, suggesting that the cells reached metabolic steady state after 5 d. In addition, the rate of glucose consumption and biomass increase was constant from 4.5 d to 5.5 d (data not shown).

Label measurements and estimation of instrument precision
Cell cultures were labelled to isotopic steady state by growth on 100% [1-13 C]glucose for 5 d under elevated O 2 and standard O 2 conditions. Quantitative 1D 13 C NMR spectroscopy was used to obtain the label measurements necessary for estimation of metabolic fluxes from soluble metabolites, protein amino acids and glucose digested from starch. The contributions to each assigned peak were analysed by line-fitting (Fig. 4, inset) and the resulting label measurements were combined appropriately to give relative cumomer abundances. This procedure yielded, from 3 biological replicates, a total of 389 relative cumomer abundance measurements for the standard O 2 condition and 429 measurements for the elevated O 2 condition. All measurements for a single biological replicate came from the same batch of cells. The complete measurement dataset is given in File S1.
To optimise flux estimates 13C-FLUX minimises the variance-weighted sum of squared differences between the experimental labelling data and simulated labelling data generated by 13C-FLUX using the flux estimates. The algorithm is therefore guided to an optimal flux solution by the size of the errors assigned to individual label measurements, and thus the accurate assignment of errors is likely to be important for the reproducible determination of the true flux solution. To accommodate this in the analysis an empirical relationship was established between relative peak error and peak signal-to-noise ratio (SNR) (Fig. 4). This formula was used to assign error estimates to individual 13 C NMR labelling measurements and hence to relative cumomer abundances.

Exploration of network structure and active reactions
The 13C-FLUX implementation of the sequential quadratic programming algorithm Donlp2 (Peter Spellucci, Technische Universität Darmstadt, Germany) was used to fit the free net and exchange fluxes to the measured isotopomer data. Fluxes were fitted to all three biological replicates for a single treatment simultaneously to give a single solution that should represent the average flux state of the three replicates. To constrain the flux solution within known bounds, we used the rates of glucose consumption and biomass accumulation, and the biomass composition, (Table 1) to calculate input and output fluxes for cell suspensions grown under elevated and standard O 2 (Sriram et al, 2006;Alonso et al, 2007a)(File S2). Fluxes derived from biomass measurements were constrained to the mean measured value and were not allowed to vary during the fitting procedure (Schwender et al, 2006). Allowing these output fluxes to vary (Alonso et al, 2007a) led to unacceptably large discrepancies between the fitted biomass fluxes and the measured values.
Free net and exchange fluxes were chosen as recommended elsewhere (Wiechert et al, 2001), and exchange fluxes were constrained to zero for steps considered to be thermodynamically irreversible. At this stage modifications to the network were incorporated to reflect the available isotopomer and biomass data, and to test whether alterations to the network structure could improve the fit of the data as measured by the sum of squared, weighted differences between the real and simulated data. The labelling data from amino acids synthesised from cytosolic (alanine, see Miyashita et al., 2007) and plastidic (valine, isoleucine and leucine) pyruvate, indicated that pyruvate in these compartments is not at isotopic equilibrium. This was modelled by introducing separate pools of pyruvate in the plastid and cytosol (ppyruvate and cpyruvate, Fig. 5). However the introduction of an irreversible uptake of pyruvate into the plastid did not improve the fit of the data suggesting that the labelling data contained little information on the uptake of pyruvate. A similar result was obtained with sunflower embryos (Alonso et al., 2007a), and a recent investigation, in which decreased plastidic pyruvate kinase produced severe decreases in lipid accumulation in Arabidopsis, is consistent with the hypothesis that uptake from the cytosol is not the main route by which the plastidic pyruvate pool is maintained (Andre et al., 2007).
Recent data also suggests that cytosolic and plastidic isoforms of NADP-malic enzyme (ME) are expressed constitutively in heterotrophic tissues of A. thaliana (Gerrard Wheeler et al., 2008).
Whilst addition of a plastidic NADP-ME (pmalic oxaloacetate ppyruvate, Fig. 5) to the model did not improve the fit of the data, analysis using EstimateStat suggested that this flux could be determined with good precision from our labelling data and it was therefore included in the final model. Addition of a cytosolic isoform of NADP-ME (oxaloacetate cpyruvate, Fig. 5) did not improve the fit above that obtained with PEPCK and the plastidic NADP-ME and in addition its presence greatly increased the flux errors predicted by EstimateStat. For this reason cytosolic NADP-ME was not included in the final model. However, in trials where cytosolic NADP-ME was included the fitting process consistently assigned it a small but significant flux, and the flux through PEPCK decreased to zero. This suggests that what we describe as PEPCK flux (the reversibility of PEPC) could equally be described as flux through cytosolic ME.

Model validation and statistical analysis
The final network structure is shown in Fig. 5, and in 13C-FLUX format, in File S3. The fitting procedure was initiated 150 times for both datasets with random initial flux estimates, and a feasible flux solution was found in over 85% of fits under both conditions. Some free fluxes converged to similar values in each run of the fitting algorithm, and for these fluxes the distribution of solutions was unimodal (Fig. 6). These fluxes included all the net fluxes, apart from flux through cytosolic aldolase under elevated O 2 conditions, and all the exchange fluxes associated with the TCA cycle. In contrast some of the remaining free fluxes did not converge to similar values in each run, with the solutions appearing to adopt a random distribution, suggesting that there is little information in the labelling data to constrain the fluxes to a particular value. The flux solution giving the lowest sum of squared weighted differences for each dataset was taken forward for further analysis and is referred to as the optimal flux solution. Figure 7 shows that there was good agreement between the observed labelling data for the standard O 2 and elevated O 2 conditions and the labelling data predicted from the optimal flux solutions. Moreover, all measurements contributing more than 1% of the total sum of squared weighted differences could be removed from the fit (24 measurements contributing 40% of the residuum for standard O 2 conditions and 21 measurements contributing 29% of the residuum for elevated O 2 conditions) without altering the optimum flux solution, demonstrating that these poorly fitting measurements were not important in constraining the fit. It can be concluded that the optimal flux solutions provide adequate descriptions of the labelling data.
The errors for the optimal fluxes were derived using EstimateStat and are summarised in Table 2.
In order to ensure that biological error present within the replicate label measurements and biomass fluxes was translated into errors in the flux estimates, the labelling data, excluding the poorly fitting data described above, were reduced to a single set of measurements (see Materials and Methods). It was then possible to use EstimateStat in combination with replicate fitting experiments ( Fig. 6) to define a list of fluxes that were statistically well determined (Wiechert et al., 2001). that elevated O 2 brings about an increase in flux through the TCA cycle and through the biosynthetic pathways associated with it. However, the proportion of carbon entering the TCA cycle that is used for biosynthesis is unaffected by the increased O 2 concentration, with 38 ± 5% of carbon used for biosynthesis under standard O 2 conditions and 33% ± 2% under elevated O 2 conditions. Moreover, if the net fluxes throughout the network are expressed relative to the rate of glucose uptake it is apparent that elevated O 2 did not bring about a major rearrangement of the metabolic network, either at the level of the TCA cycle, or at the level of the whole network ( Fig. 9).

DISCUSSION
This study aimed to quantify multiple fluxes within the central carbon metabolism of A. thaliana, and to investigate the effect of an altered O 2 concentration on the relationship between respiratory and biosynthetic fluxes around the TCA cycle. The flux maps reported here are the first to be obtained using steady state metabolic flux analysis for Arabidopsis and they show that whilst an elevated O 2 concentration in a cell culture can increase fluxes throughout the metabolic network, and alter the abundance of soluble metabolites, these changes do not require either a major reorganisation of the network or a change in the balance between respiratory and biosynthetic flux.

Reliable quantification of metabolic fluxes
Flux quantification was carried out using data obtained by 1D 13 C NMR from three biological replicates of a [1-13 C]glucose labelling experiment, leading to a set of statistically well determined flux estimates that appear to represent the global optimum flux solution. The MFA protocol incorporated several refinements aimed at improving the precision and reliability of the flux estimates.
First, to increase the likelihood of being able to accurately quantify TCA cycle fluxes, we began our investigation by predicting the optimal 13 C-labelled precursor to use for isotopic steady state labelling experiments. Whilst the approach used here considered fewer parameters than recent work in this area (Ghosh et al., 2006;Libourel et al., 2007), it nevertheless indicated that the use of [1-13 C]glucose would permit the accurate quantification of TCA cycle fluxes and associated biosynthetic and anapleurotic fluxes. In order to extend this relatively simple approach beyond the most commonly available isotopomers of glucose it will be necessary to develop methods for predicting how the precise mixture of 13 C labelled precursors influences the labelling measurements that can be made.
Secondly, in order to quantify the biosynthetic capacity of the TCA cycle more completely than hitherto, measurements of organic and amino acids derived from the TCA cycle were incorporated into the model (Fig. 2 respiration and biosynthesis. In contrast to previous work in this field all of the biomass and labelling measurements were made on cell cultures initiated from the same stock cultures and grown concurrently. This should ensure that the labelling data and biomass data are consistent with each other, which may not otherwise be the case, particularly for soluble metabolites which show significant batch to batch variability. Finally, in order to avoid the need to calculate averages of relative cumomer abundances, which are non-linear functions of metabolic fluxes, fluxes were fitted simultaneously to three separate biological replicates. This process was assisted by assigning specific errors to individual measurements on the basis of signal-to-noise ratio (Fig. 4). Biological error was applied to the estimates of flux measurement precision (Table 2)  The relative simplicity of the protocol developed in this study suggests that it may be possible to determine metabolic fluxes in heterotrophic tissue cultures of A. thaliana reasonably easily, potentially allowing a wide range of genotypes and environmental conditions to be analysed. The large datasets collected in this study will allow us, by exploiting the ability of 13C-FLUX to determine the sensitivity of individual label measurements to changes in fluxes (see Schwender et al., 2006, supplementary Fig. 6), to identify the minimum dataset needed to acquire accurate and robust estimates of metabolic fluxes in A. thaliana. Since the acquisition and analysis of large numbers of 1D 13 C NMR measurements from multiple sources (soluble metabolites, starch and protein amino acids) is currently a time-consuming process this refinement should reduce the experimental effort required for future investigations. In order to apply our approach to different genotypes it will, however, be necessary adapt the method for systems other than cell suspension cultures, which require considerable time and effort to establish.

The TCA cycle in A. thaliana
Whilst the flux distribution centred on the TCA cycle in heterotrophic A. thaliana cell suspensions under standard O 2 conditions is qualitatively similar to flux distributions in other plant species (Rontein et al., 2002;Alonso et al., 2007a, Alonso et al., 2007bSriram et al., 2007)  other tissues, the flux through NAD-ME as a percentage of pyruvate uptake was the lowest yet measured in plants (2.5%) contrasting with the previous lowest value found in tomato cell suspensions at pre-stationary phase (6%), and values as high as 66% found in oilseed rape embryos (Schwender et al., 2006). This result suggests that the dominant role of this enzyme in production of pyruvate, as proposed by others (Tronconi et al., 2008), may be a specific feature of leaves in which malate that accumulated during the light is metabolised in the subsequent dark period. Since our cells show no sign of hypoxia under standard O 2 conditions the lack of NAD-ME flux is consistent with the proposed role for NAD-malic enzyme under low O 2 (Roberts et al., 1992;Edwards et al., 1998).
In A. thaliana under standard O 2 conditions 38% of carbon entering the TCA cycle was used for biosynthetic processes, including the synthesis of protein and accumulation of amino acids and organic acids (Fig 8, Table 2). Little to no carbon was withdrawn from the TCA cycle for plastidic fatty acid synthesis, consistent with the view that the main route of carbon supply for fatty acid synthesis is via plastidic pyruvate kinase, and not via plastidic malic enzyme or uptake of pyruvate from the cytosol (Schwender et al., 2006;Alonso et al., 2007a;Andre et al., 2007). Withdrawal of TCA cycle carbon for biosynthesis has also been observed in several MFA analyses, and the proportion of carbon removed presumably reflects the balance between the demand for biosynthetic precursors, and the ATP required to convert those biosynthetic precursors into end products and maintain other cell functions. For example, in Brassica napus embryos (Schwender et al., 2006) 70% of carbon entering the TCA cycle was used for biosynthetic processes. Here, photosynthetic production of ATP reduced the need for the TCA cycle to generate reductant for the mitochondrial electron transport chain, and most of the carbon entering the TCA cycle was removed for the elongation of fatty acids in the cytosol. On the other hand, in sunflower embryos (Alonso et al., 2007a) only 6% of carbon entering the TCA cycle was used for biosynthesis, confirming the importance of the TCA cycle in meeting cellular energy demand in this system.

The metabolic response to increased O 2 concentration
The increased rates of biomass accumulation and glucose consumption in the A. thaliana cell suspension culture at elevated O 2 concentrations were associated with increases in net fluxes throughout the metabolic network (Fig. 8, Table 2). These higher fluxes corresponded to higher rates of ATP synthesis: calculations based on the optimal flux solutions, assuming that 2.5 molecules of ATP are produced for each NADH and 1.5 molecules produced for each FADH 2 (Brand, 1994), suggest that the rate of ATP production increased by around 50%, from 15 mmol d -1 flask -1 to 23 mmol d -1 flask -1 , under elevated O 2 conditions. Increased ATP levels and rates of protein synthesis have also been detected in response to increased O 2 levels in cell suspensions of peanut (Verma and Marcus, 1974). While an increase in O 2 concentration would increase ATP production if the oxygen availability under standard conditions limited the activity of cytochrome oxidase, this seems unlikely given that the O 2 levels in the cell suspension culture were well above the K m of cytochrome oxidase (0.13 µM in soybean root mitochondria (Millar et al., 1994)). Thus there would have to be a large O 2 concentration gradient between the medium and the mitochondrion for this explanation to be plausible. Moreover measurements of soluble metabolites gave no indication that the cells were hypoxic under standard conditions. For example, the abundance of alanine and GABA was significantly lower under standard conditions than under While manipulation of O 2 concentration brought about marked changes in absolute fluxes, expressing flux relative to the rate of glucose uptake showed that there was no major rearrangement of the metabolic network despite the associated changes in soluble metabolite levels (Fig. 2). In particular MFA indicates that ratios of internal fluxes in the TCA cycle and elsewhere remained almost constant (Fig. 9), a result in keeping with the unchanged biomass proportions (Table 1). Thus the changes in levels of organic acids, amino acids and sugars can be produced without major changes in relative fluxes within central metabolism. The fluxes that lead to the accumulation of soluble metabolites in this system are very small compared to fluxes though the core of carbon metabolism (Table 2), so only small net changes in flux are required to produce large changes in the relative abundance of the soluble metabolites. For example the 72% decrease in citrate accumulation relative to glucose uptake that occurred at elevated O 2 could be produced with a change in citrate synthase flux of only 4%. Similarly subtle changes in relative fluxes around OAA may be responsible for some of the changes in soluble metabolite abundances (Fig. 2). For example, under elevated O 2 conditions there was a detectable (22%) decrease in relative flux through PEPC (Figs 8 and 9) whilst flux through PEPCK increased (Table 2), possibly contributing to the decreased abundance of organic acids. However, in general, it appears that the differences in rates of accumulation of soluble metabolites arose from rearrangements of the central metabolic network that are smaller than can currently be detected using MFA.
Overall, the O 2 concentration did not exert significant control over the balance between respiration and biosynthesis under these conditions, even though it had a significant influence on the growth of the cells. Heterotrophic Arabidopsis cultures may therefore respond to changes in O 2 concentration by altering rates of respiration, biosynthesis, and ultimately growth proportionally, such that the demand for ATP to support biosynthesis is balanced by the rate of respiratory processes that generate ATP. The stability of the central metabolic network also complicates the interpretation of metabolite profiling data. The levels of most metabolites represent only a very small fraction of the total biomass making it unlikely that changes in level will reflect significant changes in flux within the central network. Thus the changes in metabolite abundance caused by altering the availability of oxygen (Fig. 2) do not lead to easily identifiable perturbations in the flux map (Fig. 8). Further improvements in the accuracy and precision of MFA may alleviate this problem, but it will also be important to complement MFA with the continued development of sophisticated models of plant metabolism (Sweetlove et al., 2008). Ultimately this can be expected to establish the relative merits of composition and flux as the basis for defining metabolic phenotypes (Ratcliffe and Shachar-Hill, 2005).

CONCLUSION
A comprehensive description of the fluxes through the TCA cycle and associated pathways in an Arabidopsis cell suspension has been obtained using a robust steady state stable isotope labelling protocol. Increasing the concentration of dissolved oxygen increased fluxes throughout the network and caused changes in the soluble metabolite profile, while at the same time having no effect on the proportion of carbon entering the TCA cycle that was used for biosynthesis, and no significant impact on the relative flux distribution in the central network. Ultimately, while the mechanism by which the cells respond to increased oxygen has yet to be established, the study demonstrates the utility of MFA as a tool for probing the impact of an environmental perturbation on the operation of the central metabolic network.

Oxygen electrode measurements
Measurements of dissolved oxygen concentration were made using a Clark type oxygen electrode at 21 o C. Arabidopsis cell suspension (1 ml) was transferred to the electrode chamber and oxygen consumption was monitored until the trace became linear. The linear portion was extrapolated to time zero to determine the oxygen concentration at the point of sample addition.

Isotopic steady state labelling
Cells were labelled to isotopic steady state by subculturing light grown cell suspension into medium where a proportion of the unlabelled glucose was replaced with 13 C-labelled glucose (Cambridge Isotope Laboratories, Andover, Massachusetts and Sigma-Aldrich, Gillingham, U.K.). This procedure has been shown to have no discernible effect on the flux distribution through the Arabidopsis metabolic network (Kruger et al., 2007b). Biological variation was assessed by subculturing from 3 separate light grown stocks that had been maintained independently for several weeks. Cultures were incubated on an orbital shaker in the dark at 21 o C for 4.5, 5 or 5.5 d as appropriate. Cells were harvested by vacuum filtration through a single paper filter, washed with 210 ml glucose-free growth medium, weighed, and immediately frozen in liquid N 2 . Tissue was stored at -80 o C prior to analysis.
Soluble metabolites were extracted from frozen tissue labelled to isotopic steady state using perchloric acid (Kruger et al., 2007b). Following the final freeze drying step samples were redissolved in 10% 2 H 2 O, with 10 mM EDTA, 25 mM 1,4-dioxane and 10 mM KH 2 PO 4 /K 2 HPO 4 pH 7.5 for NMR spectroscopy. Starch was extracted from the insoluble residue remaining from perchloric acid extractions. The residue was washed, and autoclaved in 100 mM sodium acetate (pH 4.8) for 2 h. Gelatinized starch was then enzymatically digested overnight at 37 o C with 30 units of α -amylase (Roche, Basel, Switzerland) and 5 units of amyloglucosidase (Roche). The supernatant, containing glucose released from starch, was freeze dried and redissolved in 10% 2 H 2 O with 25 mM 1,4-dioxane for 13 C NMR analysis.
Protein was extracted by repeated washing of ground, lyophilized tissue with phosphate buffered saline (130 mM NaCl, 100 mM Na 2 HPO 4 /NaH 2 PO 4 pH 7.0). Prior to hydrolysis, protein was precipitated using 12% TCA, washed with ice-cold acetone, and resuspended in 6 M HCl.
Hydrolysis was carried out in Pierce hydrolysis tubes (Pierce, Rockford, Illinois): samples were degassed and flushed with N 2 three times then heated at 95 o C for 24 h under vacuum. Samples were freeze dried to remove HCl, and redissolved in 10% 2 H 2 O with 25 mM 1,4-dioxane, pH 7.5 for 13 C NMR analysis.

Biomass analysis
Growth rate of cell suspensions was determined from fresh weight recorded during harvest, and converted to change in dry weight by assuming that fresh cells contained 95% water by weight. This value was supported by experiments in which cell mass was determined before and after freeze drying. Measurements of the abundance of protein, amino acids, cell wall, starch, lipids and soluble metabolites were made using either tissue labelled to isotopic steady state with [1-13 C]glucose or tissue grown concurrently from the same stock cultures.
Protein extracted with phosphate buffered saline was quantified using the Bradford assay. The amino acid content of labelled protein hydrolysates was determined by HPLC (Bruckner et al., 1995)  separated using a reverse phase C18 column and quantified by fluorescence using standard curves.
Starch was quantified by autoclaving duplicate samples of ground, unlabelled, lyophilized tissue for 1 h in 25 mM sodium acetate pH 4.8. Fifteen units of α -amylase (Roche) and 5 units of amyloglucosidase (Sigma-Aldrich) were added to one of the samples and both samples were placed at 37 o C overnight. Glucose in the supernatants from both samples was quantified using a spectrophotometric assay (Sweetlove et al., 1996) and the abundance of starch was calculated from the difference in the amount of glucose in the enzyme treated and untreated samples. The same spectrophotometric assay was used to quantify glucose in the growth medium.
Cell wall was extracted by repeated washing of a known mass of unlabelled ground lyophilized tissue with a mixture of phenol, acetic acid and water in the ratio 2:1:2 (Sriram et al., 2006).
Insoluble material remaining was washed with distilled water to remove residual phenol, freeze dried and weighed. The mass was corrected for the presence of contaminating starch using the starch quantification protocol above.
Lipids were extracted from a known mass of ground, labelled, lyophilized tissue using hexane and isopropanol according to an established protocol (Hara and Radin, 1978;Mhaske et al., 2005).

Solvent was removed by gentle heating and lipids quantified by weight.
Soluble metabolites were extracted with methanol according as described elsewhere (Le Gall et al., 2003). A known mass of ground, lyophilized tissue was shaken with 70 % methanol, 30 % buffer (1 mM sodium 3-trimethylsilylpropionate 2 H 4 (TSP), 1 mM EDTA, 50 mM KH 2 PO 4 , 50 mM KH 2 PO 4 ) for 30 min. Solvent was removed under vacuum and samples were redissolved in 800 µl 2 H 2 O for 1 H NMR analysis. Quantification was based on the addition of a known quantity of TSP to representative samples.

NMR spectroscopy
Spectra were recorded on a Varian Unity Inova 600 spectrometer (Varian Inc, Palo Alto, California). 1D 13 C NMR spectra were recorded at 150.9 MHz using either a 10 mm broadband or a 10 mm 13 C/ 31 P switchable probe, and in all experiments Waltz16 decoupling was applied during the detection period to decouple 1 H signals. Spectra were referenced to 1,4-dioxane at 67.3 ppm.
Where absolute quantification of labelling was not required a recycle delay of 6 s was used, and the NOE was induced during the relaxation delay to increase SNR values. For quantitative data the recycle delay was extended to 19 seconds and the NOE was not induced. Acquisition times of the order of 60 h (10,240 scans) were required to obtain suitable spectra from labelled samples for accurate line-fitting. Spectra were acquired in blocks of 1,024 scans that were manually summed following inspection to confirm that there was no significant degradation during acquisition. 1D 1 H NMR spectra were recorded at 600 MHz using a 5 mm HCN triple resonance probe and the standard Varian pulse program. Presaturation was applied during the relaxation delay to suppress the water signal and spectra were referenced to TSP at 0 ppm. 2D 1 H/ 13 C gHMBC and gHSQC spectra were recorded at 600 MHz ( 1 H) and 150.9 MHz ( 13 C) using a 5 mm HCN triple resonance probe and standard Varian pulse programs. WURST-40 or Garp1 decoupling was applied during the detection period to remove 13 C-1 H coupling.
All spectra were processed and analysed using NUTS (Acorn NMR). 1D 13 C spectra were processed using a line-broadening of 2.5 Hz for spectra requiring line-fitting or 1 Hz for spectra requiring integration. 1D 1 H spectra were processed using a line-broadening of 1 Hz. 1 H and 13 C assignments were based on literature values, comparison with pure standards and the results of 2D NMR experiments. Spectral deconvolution (line-fitting) of 1D 13 C spectra was carried out using the line-fitting subroutine in NUTS. During line-fitting, resonance frequency, signal intensity, linewidth and fraction Lorentzian lineshape were varied to minimise the difference between the real and simulated spectra.
The SNR values for the 1D 13 C NMR signals from soluble extracts of cells labelled to isotopic steady state with [1-13 C]glucose varied over more than two orders of magnitude between different metabolites and different experimental conditions. Variation in extraction efficiency and sample fresh weights also contributed to considerable variation in signal intensity for the same metabolites between biological replicates. Since instrumental precision depends on SNR, it would have been incorrect to assign the same relative error to every label measurement during optimisation of fluxes. We therefore recorded replicate quantitative 1D 13 C NMR spectra of standard samples of organic and amino acids (25 mM alanine, 50 mM citrate, 5 mM glutamate, 0.5 mM aspartate and 1 mM malate) and [U-13 C 6 ]glucose. Signal intensities were measured by line-fitting and the relative error of the same signal over three replicate spectra was correlated with the corresponding SNR ( Fig. 4). An empirical relationship between SNR and relative peak error was determined (Fig. 4) and this formula was used to assign error estimates to 13 C NMR labelling measurements.

Metabolic modelling
Metabolic modelling was carried out using 13C-FLUX (version 20050329) (Wiechert et al., 2001 and references therein). The EstimateStat component of 13C-FLUX was used to refine the metabolic network. To do this initial flux estimates were taken from the literature (Rontein et al., 2002) and lists of 1D 13 C NMR measurements were taken from steady state experiments carried out previously in our lab (P. Lelay, unpublished observations). The CumoNet component of 13C-FLUX was then used to predict a set of cumomer and isotopomer abundances consistent with the initial flux estimates, and these abundances were combined appropriately to produce pseudo 1D 13 C NMR datasets in which the label measurements perfectly reflected the underlying fluxes. An error of 5% was assumed for each measurement within this dataset. EstimateStat was then used to determine the errors associated with the flux estimates. The same approach was used to predict the optimum labelled precursor for minimisation of predicted flux estimate errors. Here the refined network was used in conjunction with sets of 1D 13 C NMR measurements derived from experiments specific to the fed precursor.
The estimated flux errors produced by EstimateStat are sensitive to measurement configuration and relative measurement errors, but do not depend on absolute label measurements (Wiechert et al., 2001). To derive flux error estimates that incorporate biological error we determined the relative error in cumomer abundances between the different biological replicates. To do this we first normalised the labelling data from the three biological replicates using the group scaling factors estimated during the fitting process (Wiechert et al., 2001). The average and standard deviation of replicate normalised cumomer abundances were then calculated, and these measurements used with EstimateStat to estimate flux errors. In order to combine the errors in fluxes derived from the biomass data with errors associated with the labelling data, the biomass derived fluxes were set as free fluxes for this part of the analysis, and the biomass derived fluxes and their standard deviations were defined in the Flux Measurements section of the input file.

Statistical analysis
All indications of statistical significance are based on a Student's t-test with P < 0.05 unless otherwise indicated. The effect of labelled precursor on predicted relative flux errors (flux error/flux). Simulations used the network illustrated in Fig. S1. A relative error of less than one (indicated by the dashed line)

Supplemental
indicates that the predicted flux error is less than the magnitude of the flux itself. Abbreviations: CO 2 (CO 2 output), cALD (cytosolic aldolase), pALD (plastidic aldolase), PK (pyruvate kinase), G6PDH (glucose 6-phosphate dehydrogenase), PDH (pyruvate dehydrogenase), SDH (succinate dehydrogenase), MDH (malate dehydrogenase), PEPC (PEP carboxylase), ME (mitochondrial NAD-malic enzyme).   The dependence of relative error of peak areas on peak SNR in 1D 13 C NMR spectra. From these data power regression was used to define an empirical relationship between relative error (RE) and SNR: RE = 0.62(SNR) -0.66 . Relative error was calculated as the standard deviation of three www.plantphysiol.org on August 14, 2017 -Published by Downloaded from Copyright © 2008 American Society of Plant Biologists. All rights reserved. measurements made on the same sample divided by the average peak area. Measurements were made on a sample of [U-13 C 6 ]glucose with a low SNR and on a mixture of alanine, citrate, glutamate, aspartate and malate at higher SNR. Errors are ± one standard deviation (n = 3). Inset, line-fitting of a representative peak from a 1D 13 C NMR spectrum. Line-fitted peaks are indicated with dashed lines. Peak annotations indicate the cumomer abundance defined by the area of the line-fitted peak; 1 corresponds to the presence of 13 C at positions 1 to 6 in citrate (standard numbering), 0 corresponds to the presence of 12 C and x corresponds to either 13 C or 12 C.     Table 2  Reactions that carry no flux under the specified condition are not included. Forward and reverse fluxes for ana1 (PEPCK) are illustrated. Exchange fluxes for tca4 and tca5a+5b were also determined ( Table 2) but are omitted for clarity. The percentage of carbon entering the TCA cycle used for biosynthesis was calculated as: (4ana3+4asp+5glu+6scit+4smal+4ssucc+4sfumOUT)/(3cmex+4ana1) Figure 9.
The effect of elevated O 2 on metabolic fluxes relative to the rate of glucose uptake. Fluxes corresponding to the optimal solution are divided by the rate of glucose uptake, thereby permitting changes in organisation of the metabolic network to be distinguished. Errors are ± one standard deviation as determined using EstimateStat. Flux names correspond to those given in Fig. 5 and    Col 5 Col 7 Col 9 Figure 1. The effect of labelled precursor on predicted relative flux errors (flux error/flux). A relative error of less than one (indicated by the dashed line) indicates that the predicted flux error is less than the magnitude of the flux itself. Simulations used the network illustrated in Fig. S1. Abbreviations: CO 2 (CO 2 output), cALD (cytosolic aldolase), pALD (plastidic aldolase), PK (pyruvate kinase), G6PDH (glucose 6-phosphate dehydrogenase), PDH (pyruvate dehydrogenase), SDH (succinate dehydrogenase), MDH (malate dehydrogenase), PEPC (PEP carboxylase), ME (mitochondrial NAD-malic enzyme).
Relative Flux Error [1][2][3][4][5][6][7][8][9][10][11][12][13]   . Errors are ± one standard deviation (n = 2 or 3). "a" indicates statistically significant differences (Student's Ttest, p < 0.05) in metabolite abundances between the two conditions for that time point. "b" indicates statistically significant differences (Student's T-test, p < 0.05) in metabolite abundance under a particular condition between the indicated time point and 4.5 d growth. No significant differences in abundance were detected between 5 and 5.5 d growth.   Dashed boxes indicated where the subcellular localisation of a metabolite or reaction cannot be inferred from the data or from the literature. The letters "p", "c" and "m" preceding metabolite names indicate separate pools of that metabolite in the plastid, cytosol and mitochondrion respectively. PPP indicates the pentose phosphate pathway. Output of CO 2 from the system is included in the model but not illustrated here.    Table 2 and the number of carbon atoms in the metabolites involved in each reaction. Reactions that carry no flux under the specified condition are not included. Forward and reverse fluxes for ana1 (PEPCK) are illustrated. Exchange fluxes for tca4 and tca5a+5b were also determined (   Figure 9. The effect of elevated O 2 on metabolic fluxes relative to the rate of glucose uptake. Fluxes corresponding to the optimal solution are divided by the rate of glucose uptake, thereby permitting changes in organisation of the metabolic network to be distinguished. Errors are ± one standard deviation as determined using EstimateStat. Flux names correspond to those given in Fig.  5 and Table 2. Flux relative to glucose uptake (arb units)