Abstract

Engineering plant primary metabolism is currently recognized as a major approach to gain improved productivity. Most of the current efforts in plant metabolic engineering focused on either individual enzymes or a few enzymes in a particular pathway without fully consider the potential interactions between metabolisms. More and more evidences suggested that engineering a particular pathway without consideration of the interacting pathways only generated limited success. Therefore, a long term goal of metabolic engineering is being able to engineer metabolism with consideration of the effects of external or internal perturbation on the whole plant primary metabolism. In this paper, we developed a constraint-based model of C3 Plant Primary Metabolism (C3PMM), which is generic for C3 plants such as rice, Arabidopsis, and soybean. The C3PMM was first combined with transcriptome data to demonstrate that there is substantial coordination of mesophyll primary metabolism at both transcriptome and metabolism levels in response to elevated CO2 concentration. Secondly, maximizing CO2 uptake is a plausible target function for metabolism of a typical C3 mesophyll cell. Finally, C3PMM predicted a decrease in nitrate assimilation flux coordinated with photorespiration, and increase in ammonium assimilation flux, when CO2 concentration increases. This suggested a potential mechanistic linkage between differences in the response of ecosystems differing in nitrogen source to elevated CO2. In conclusion, all the coordination of C3 plant primary metabolism at both transcriptome level and metabolism level, and the coordination between nitrate assimilation and photorespiration under elevated CO2 concentration, are beneficial for maximized CO2 uptake rate.

Keywords

Introduction

Plant metabolic network is inherently complex. Different metabolic pathways in the network are connected with each other and influence each other. For example, in C3 plants, the nitrogen metabolism is linked to the photorespiration pathway, as demonstrated by the decrease of nitrogen assimilation under elevated CO2 [1]. Similarly, many intermediates in the Calvin cycle are shared with the pentose phosphate pathway. In the C4 photosynthetic metabolism, the four carbon acid shuttle has been hypothesized to be linked to the nitrogen reassimilation in the senescence leaf [2]. All these examples demonstrated that the photosynthetic process is integrated with many other plant primary metabolic pathways. This complex interaction between pathways determines the actual performance of the system, such as responses under changed environmental conditions. This complexity also renders the control over the fluxes of the system shared by many enzymes in the system, and as a result, makes it rather difficult to predict the results of engineering a metabolic system through single gene approach. This is reflected by the limited success of engineering for improved photosynthetic energy conversion efficiency, though various approaches to manipulate photosynthesis has been proposed [3,4]. Similarly, when chloroplast NADP-Dependent Malic Enzyme (NADPME) was transferred into rice, aberrant agranular chloroplast without thylakoid stacking was generated, with decreased level of chlorophyll content and photosystem II activity [5]; these changes might be related to potential influence of manipulation of NADP-ME on the redox state of the chloroplast [5].

Being able to engineering metabolism with consideration of the interaction between different components, either being enzymes or pathways, hold great potential to expedite the progress of metabolic engineering. To realize this, a number of fundamental questions however have to be addressed. Firstly, are the responses of metabolic pathways to internal or external perturbations coordinated? In Escherichia coli, metabolism responds to external perturbations in a highly coordinated manner [6]. However, whether there is coordination between different pathways of the plant primary metabolism has not been systematically evaluated so far. Secondly, if the pathways in plant primary metabolism respond to perturbations in a coordinated manner, do they respond in a manner to gain a particular biological function? In Escherichia coli, by comparing the flux distribution calculated using genome-scale model with the gene expression information from microarray, Edwards et al. demonstrated that the gene expressions respond to different environments in a coordinated manner to ensure maximal biomass production [6]. With plant being a multi-cellular eukaryotic system, will maximizing biomass production still be its default objective function? There have been some studies on plant gene expression response to environment changes. For example, Liu et al. analyzed gene expression changes associated with boron deficiency in citrus to provide further information for understanding the mechanism of the different responses of scion-rootstock combinations to boron-deficiency stress [7]. D'Esposito et al. further combined genome, transcriptome, metabolome and sensorial data of three tomato varieties to quantify the transcriptional response to environmental cues, and measure the metabolic activity, demonstrated the principal role of the cell wall in fruit quality [8]. But few studies deal with the coordination of different pathways within the framework of complex plant metabolic network.

Constraint-based metabolic model has been used as an effective tool to study the phenotypic responses of organism under varying conditions. It is a formal and mathematical representation of reconstructed metabolism as a genome-scale network, which consists of collections of metabolic reactions, their stoichiometry, the enzymes and the genes that encode them [9]. In contrast to kinetic models of metabolism [10-12], constraint-based models can predict steady-state fluxes in large-scale metabolic network by using Flux Balance Analysis (FBA) without requiring detailed kinetic parameters. The principle of FBA is based on stoichiometric matrix, by adding some thermodynamic information, and constraints of metabolism capacity. FBA identifies the optimal state of the metabolic network that would allow the system to achieve a particular objective, typically the maximization of production of biomass or specific target compound in an organism [13]. So far, constraint-based models have been developed for many prokaryotic systems [14], and even eukaryotic systems, such as Saccharomyces cerevisiae [15], mouse [16,17] and human [18,19]. Though there are many inherent difficulties to develop constraint-based models of plant metabolic network because of high compartmentalization, numerous reversible reactions, and unknown biochemical steps [20-22], some constraint-based models have already been developed in plants such as Arabidopsis [23-25], barley [26], rice [27] and C4GEM (a generic model for C4 plants maize, sorghum and sugarcane) [28]. In our previous work, we have refined AraGEM and C4GEM models to compare the effects of gene deletions on photosynthesis and biomass [29]. However, there are still challenges in refining the genome-scale plant model, even for the more studied Arabidopsis, such as filling the gaps, adding transport reactions, modifying particular constraints, and updating the Gene- Proten-Reaction rules. To objectively answer the question whether plants respond to perturbations in a coordinated manner, we chose to focus on the relative clear primary metabolism of C3 plant. It has been reported in Yeast that there is correlation between gene expression and metabolic flux in some primary pathways, while no correlation in secondary metabolism [30]. Therefore, it is more advisable to evaluate the coordination between the two levels using the primary metabolism model of C3 plant. Furthermore, since the biochemical reactions in primary metabolism of most C3 plants are conserved, such as rice, Arabidopsis, and soybean, we can apply the generic model to address questions common to different C3 plants.

In this study, we aim to uncover the coordination of plant primary metabolism in response to climate change. Firstly, we developed a constraint-based model of the plant primary metabolic network in C3 mesophyll cell (C3PMM); Secondly, tested whether there is coordination between flux change and gene expression change under elevated/ambient CO2 concentrations; Then, we used the C3PMM to study the interaction between three major functional modules in the mesophyll cell: The Calvin cycle, photorespiration pathway, and nitrogen metabolism.

Methods

Reconstruction of the constraint-based model of primary metabolism in a C3 mesophyll cell The C3 Plant Primary Metabolic Model (C3PMM) is compartmentalized into cytosol, mitochondria, chloroplast and peroxisome. Initial pathway and reaction information were obtained from the databases including KEGG, PlantCyc and Brenda. For the compartmentalization, we assigned the basic reactions according to biochemistry knowledge, for example, the Calvin cycle in chloroplast, TCA in mitochondria. Also some reactions’ location is determined according to KEGG and PlantCyc. For those reactions with no certain annotation, we manually curated based on the relative biochemical and physiological literatures. We reconstructed and modified the model using a visualization toolbox SimBiology in MATLAB. The reconstructed model includes photosynthesis, photorespiration, Glycolysis, TCA cycle, amino acid metabolism, starch synthesis, sucrose synthesis, and other important transport processes, totally 134 reactions and 131 metabolites, as shown in Table 1. The detailed reaction list (Tables S1-S3) and SBML model file (Model S1) could be downloaded from the supplemental data.

The Flux Balance Analysis is a constraint-based metabolic network analysis method, by adding some stoichiometric information, thermodynamic information, and constraints of enzyme capacity. The behavior of the system can be restricted in a closed solution space, and then the optimal solution can be calculated by linear programming. The mass conservation and steady state assumption is necessary:

Where S is the stoichiometric matrix of metabolic system and si, j is the stoichiometric coefficient of Metabolite i in reaction j, v j is the flux of each reaction in the metabolic system. In addition to mass balance, a set of constraints are required to restrict the feasible solution space, which can be defined as below:

Then the optimal metabolic flux distribution can be found using linear programming given an objective function as following:

The flux variability analysis can calculate the full range of numerical values for each reaction flux in a network. This is carried out by optimizing for a particular objective, while still satisfying the given constraints in the system.

The different light intensities were simulated by setting flux value of the input photon, which is the reaction named ‘Ex-photon’.

3Pi + 3ADP + 2H+ + 2NADP + 8 photon⇒3ATP + 2NADPH

Then the flux of light-reaction will be 1/8 of the flux of ‘Ex-photon’, which can affect the flux distribution of other reactions including CO2 fixation by Rubisco. By this means, the optimal CO2 assimilation rate can be predicted at different light intensities.

Another method to characterize the solution space of metabolic flux is uniform random sampling, which does not require the objective function. Specifically in sampling, we used an implementation of the Artificial Centered Hit-and-Run (ACHR) sampler algorithm with slight modification to generate such a set of flux distributions that uniformly sample the space of all feasible fluxes [31,32]. Initially, a set of 5000 nonuniform pseudo-random points, called warm-up points, was generated. In a series of iterations, each point was randomly moved while keeping it within the feasible flux space. This was accomplished by choosing a random direction, computing the limits on how far a point could travel in that direction (positive or negative), and then choosing a new point randomly along that line. After numerous iterations, the set of points was mixed and approached a uniform sample of the solution space and 2000 points was loaded for analysis. The set of flux distributions obtained from sampling can be interrogated further to answer a number of questions related to the metabolic network function. Here, we used the results of sampling to predict the correlated reaction sets, where the reactions in each set show similar changes under perturbation. Also we made 2000 samplings and used the maximal flux value of CO2 assimilation by Rubisco at each CO2 concentration to build the A-Ci curve. The FBA, FVA, and sampling analysis of our model were all conducted using COBRA toolbox.

Description of the elevated/ambient CO2 conditions

Because Rubisco catalyzes two reactions with CO2 and O2 respectively, it is not suitable to simulate the optimal CO2 fixation rate without any constraints on oxygenation. In our C3PMM, we combined the two reactions catalyzed by Rubisco into one, and set the ratios of rate between RuBP carboxylation and oxygenation, as following:

Where RuBP is D-ribulose 1,5-bisphosphate, 3PGA is 3-phospho- D-glycerate, PGCA is Phosphoglycolate, Rubico is ribulosebisphosphate carboxylase.

The r value represents the ratio of carboxylation and oxygenation, which can be calculated according to the well-known Farquhar’s model. The different CO2 concentrations were represented by setting different ratios (r value) of the rate between CO2 (Vc) and O2 (VO) (Table 2) [29].

CO2 concentration (ppm)

Vc/Vo (r)

380

4.33

550

6.26

800

9.11

1000

11.39

Table 2: The ratio of carboxylation (Vc) and oxygenation by Rubisco under different CO2 concentrations.

Gene expression datasets

The microarray dataset for 800ppm CO2 was retrieved from TAIR database (Accession=HybData: 1005823716), which measured the gene expression of Arabidopsis thaliana Columbia-0 grown at 800ppm CO2 for 30 days. The positive fold change is the ratio of expression signal at 800 vs 380ppm CO2. The genes with ratio higher than 1.5 were labeled as up-regulated (according to the description of the microarray dataset). Similarly, the negative fold change represented decreased expression, which is the ratio between 800 vs. 380 ppm CO2, so we set genes with fold change less than -1.5 as down-regulated. The microarray dataset for 1000 ppm CO2 was also retrieved from TAIR database (Accession=HybData: 1005824007), we determined the upand down-regulated genes in the same way as the dataset for 800 ppm CO2. The other microarray dataset for 550 ppm CO2 was obtained from Li et al. We selected the average expression ratio of Arabidopsis thaliana Columbia-0 (Col-0) and Cape Verde Island (Cvi-0) harvested on June 21, 2005. The significantly up and down regulated genes between 550 and 380 ppm CO2 are determined using the cutoff fold change of 1.3 and 0.77 respectively as mentioned in the reference.

Results

The constraint-based model of C3 plant primary metabolism

Plant metabolism features high degree of compartmentalization, there are multiple pools of metabolites, and the existence of duplicate pathways in different organelles. Following the protocol of constraintbased model construction [31]; we made a de novo reconstruction of C3 primary metabolic network with 134 reactions and 131 metabolites. This model includes four cellular compartments, Chloroplast [Ch], mitochondria [M], peroxisome [P], and cytosol [C], Table 1 showed the composition of the network. The flux unit in our model is μmol. kgDW-1.min-1. Since it is a generic primary metabolism model for C3 plants, the particular gene information of particular species were not specified in the Systems Biology Markup Language (SBML) file of the model. The diagram of whole network C3PMM (Figure S1) and the SBML description of the model (Model S1) are available from the supplemental data.

The topology structure of the primary metabolic network was illustrated in Figure 1, with the node size proportional to the degree of each reaction node. There are 50 hub nodes with degree larger than 10, among which 52% are predicted to be essential reactions for photosynthesis, which means when deleting these reactions the flux of CO2 fixation will be zero. In addition, we used uniform random sampling to identify the correlated sets with all reactions having correlated flux with each other across 2000 samples [32], see detail in Methods section. We found the largest correlated reaction set (green nodes in Figure 1) was significantly enriched in essential reactions (Hypergeometric test p-value=6.2e-10), with all 28 reactions being essential for photosynthesis, which indicated this part is dominant to determine the response of C3 primary metabolism to genetic or environment changes (Among the other 28 essential reactions, only 4 reactions are correlated).

Figure 1: The topology structure of C3 plant primary metabolic network. Each node represents a reaction, if one product of reaction1 is the reactant of reaction 2, there will be an edge between the two nodes. The node size is proportional to the degree of each reaction, and the label on each node represents it is included in the corresponding correlated reaction set. The labels with reaction names represent those reactions not belonging to any correlated set. The red and green nodes are essential reactions for photosynthesis, and the green nodes are also involved in the largest correlated reaction set. The blue nodes are not essential reactions for photosynthesis.

Model validation with physiological data

The curve of CO2 assimilation versus CO2 concentration (usually called A-Ci curve) is a basic measurement for model validation. In a typical photosynthetic physiological experiment, the A-Ci curve is taken usually using a saturating light level. Under this condition, the CO2 assimilation rate is indeed the maximal CO2 uptake rate at a particular CO2 concentration. Therefore, in order to match the simulation with the typically measured A-Ci curve, we used the uniform random sampling method and took the maximal CO2 uptake rate among 2000 samples under different CO2 concentrations. The predicted A-Ci curve mimics the experimentally measured A-Ci curve [33,34], as shown in Figure 2A. In order to compare the flux unit (μmol.kg-1.min-1) in our C3PMM simulation with the unit of CO2 assimilation rate by experiment (μmol.m-2.s-1), the Specific Leaf Area [35] is required. Since SLA varies with different CO2 concentration and nitrogen supply, for example, 60-150 cm2.g-1 for Arabidopsis [36], conversion from flux to real rate should be divided by a number between 4-10. Therefore in Figure 2, the predicted flux of CO2 fixation is comparable with the experimental CO2 assimilation rate. Furthermore, we predicted the optimal CO2 assimilation rate under different light intensity. The predicted gradual increase of assimilation rate with increase of light levels is also consistent with experimental observation [34], as shown in Figure 2B. In addition, we predicted response curves of CO2 uptake rate to light under different CO2 levels (Figure S2). All these predictions are consistent with experimentally observed responses. Furthermore, we predicted the gradual increase in the flux of biomass production under different CO2 levels as well, with coefficients for 10 biomass components referencing the study by Poolman et al. (Figure 3) [23]. All these results suggested that the model can reflect basic features of primary metabolism of C3 plants.

Figure 2: The validation of predicted flux of photosynthetic CO2 fixation vs intercellular CO2 concentration and light intensity. A. Relationship between photosynthetic CO2 assimilation vs intercellular CO2 concentration by uniform random sampling analysis (A-Ci curve). Under each intercellular CO2 concentration, the maximal rate of CO2 assimilation was selected from 2000 sampling results. The line with triangle represents the predicted CO2 assimilation flux by C3PMM, the dotted line is the experiment results of CO2 assimilation rate [30]. B. The predicted responses of photosynthetic CO2 uptake rate to light intensity. The line with triangle is the change of optimal CO2 uptake rate with light intensity predicted by C3PMM; the dotted line represents experimental results [30]. The unit of predicted flux is μmol.kg-1.min-1, the unit of CO2 assimilation rate by experiment is μmol.m-2.s-1, so the conversion is 60/SLA. For example, SLA in Arabidopsis is 60-150 cm2.g-1 [32], conversion from flux to real rate is divided by a number between 4-10.

Figure 3: The predicted flux of biomass production to CO2 by FBA. The optimal rate of biomass synthesis is computed under each intercellular CO2 concentration

Coordination between changes in metabolic fluxes and gene expression under ambient/elevated CO2 concentrations

In this study, we simulated different CO2 concentrations by setting different values of ratio between carboxylation and oxygenation catalyzed by ribulose-bisphosphate carboxylase (Rubisco), see Methods section. We tested five potential objective functions, including maximization of CO2 uptake, maximization of starch synthesis, maximization of sucrose synthesis, maximization of biomass production and minimization of fluxes through Glycolysis. We simulated the flux distribution using flux balance analysis under three elevated CO2 concentrations, 550/800/1000 ppm, and compared the results with the ambient CO2 concentration of 380 ppm. Results suggested that flux of photorespiration decreased more rapidly than the increase of CO2 fixation flux with increasing CO2 concentrations. Same flux distribution patterns were observed with either maximizing sucrose or minimizing Glycolysis as the objective functions.

We compared the changes in flux distribution under two elevated CO2 concentrations (800 and 550 ppm) with the corresponding gene expression changes for Arabidopsis, using the values at 380 ppm as control. The significantly up and down regulated genes between elevated and ambient CO2 were determined using the cutoff fold change by the original studies of the microarray datasets (Methods) [37]. Similar with the up- and down-regulation for gene expression, we set the reactions with flux ratio more than one fold standard deviation away the mean as differentially up- and down-regulated. Since there are different isoforms for some enzymes, we selected the genes with majority of them being up- or down-regulated, and then pick the highest absolute value of expression to represent the particular enzyme. The heatmap of predicted flux ratios and gene expression ratios between 550 and 800 ppm CO2 and ambient condition was illustrated in Figure 4. Results showed that most of the predicted fluxes change in the same trend with the expression levels (Table 3). This is especially the case when the objective function was set as maximization of photosynthetic CO2 uptake rate, there are 81.67% and 88.41% enzymes having same changes on flux and gene expression respectively at 800 and 550 ppm CO2. To ensure such correspondence between gene expression and metabolic flux is not an artifact by the particular FBA solution, the flux variability analysis was also conducted to get the range of reaction fluxes. We found most of the flux variations also change in the same trend as the corresponding gene expression level, especially there are 81.8% and 97.1% enzymes having same changes on flux range and gene expression respectively at 800 and 550 ppm CO2, when the objective function is CO2 uptake rate (Table S1).

Objective functions

Percentage of reactions with Rf consistent with Re (800ppm CO2)

Percentage of reactions with Rf consistent with Re (550 ppm CO2)

Max CO2 fixation rate

81.67%

88.41%

Max starch synthesis

76.67%

85.51%

Max sucrose synthesis

71.67%

82.61%

Glycolysis

71.67%

82.61%

Max Biomass production

46.67%

56.52%

Rf represents the ratio of flux between elevated (either 800 or 550 ppm) and ambient CO2 concentrations; Re represents the ratio of gene expression between elevated (either 800 or 550 ppm) and ambient CO2 concentrations. Reactions with Rf and Re being both up-regulated and down-regulated simultaneously are labeled as consistent.

Table 3: The consistencies of changes in the gene expression and the fluxes under different CO2 concentrations predicted with different objective functions

Figure 4: The heatmap illustration of log2 ratio in reaction flux and gene expression at 550 and 800 ppm CO2 over ambient condition. The range of the bar is different in the two panels because the distribution of gene expression values is different in the two microarray datasets.

Interestingly, all of the differentially expressed genes (under elevated CO2) are critical for maintaining the flux through the Calvin cycle. When these genes are assigned to have a flux of 0, the total CO2 uptake rates will become 0. Furthermore, under elevated CO2, the expression level of five enzymes were not identified in the microarray experiments (Arabidopsis, 800 vs. 380 ppm CO2), including 6-phosphofructokinase, diphosphate-fructose-6-phosphate 1-phosphotransferase, fructokinase, UTP-glucose-1-phosphate uridylyltransferase, and pyruvate decarboxylase. Interestingly, the reactions catalyzed by these five enzymes also showed zero fluxes under both control (380 ppm) and elevated (800 ppm) CO2 concentrations.

We found that reactions in the same pathway demonstrated similar magnitude of changes in the flux under elevated versus ambient CO2 concentrations, when the objective function was maximization of CO2 uptake rate. Figure 5 showed the ratios between fluxes of the same reaction in TCA cycle, Calvin cycle and starch synthesis, and glycolysis and sucrose synthesis. The flux changes in all pathways were provided in supplemental Table S2. For example, at 800 ppm CO2, the fluxes in photorespiration and nitrate assimilation significantly decrease (0.519, Table S2); the fluxes in light-reactions and TCA cycle show a little decrease (Figure 5A and Table S2); the fluxes in Glycolysis and sucrose synthesis show no change (Figure 5C). Furthermore, the reactions in the Calvin cycle were divided into two different stages, with one stage covering the reactions from RuBP carboxylation till the site of the fructose 6 phosphate formation; while the other stage including reactions from formation of fructose 6 phosphate to the formation of Ribulose Bisphosphate (RuBP). Under elevated 800 ppm CO2, flux in RuBP increase slightly with a ratio of 1.019 (Figure 5B), and reactions associated with the regeneration of RuBP are predicted to be slightly decreased (Figure 5B), which might reflect a decreased demand for RuBP under elevated CO2 conditions. We also found the increased flux in Aspartate and Alanine metabolism (1.55), the decreased flux in pyruvate metabolism (0.96), and no change in ethanol metabolism and hextose shuttle. Furthermore, the shift in the range of fluxes by FVA is also similar within the particular pathway, as shown in Table S1. Therefore, the primary metabolism of C3 plants responses to environmental changes in a coordinated way.

A) The coordinated fluxes in Calvin cycle and starch synthesis in chloroplast. B) The coordinated fluxes in TCA Cycle in mitochondria. C) The coordinated fluxes in Glycolysis and sucrose synthesis in cytosol. Each block represents an enzyme with EC number, each circle represents a metabolite, and arrows show the reaction flow between metabolites. The blocks with numbers on right side represent the flux ratio of corresponding pathway with same color under 800 vs. 380 ppm CO2

Flux of nitrate assimilation depends on photorespiration and decreased with elevated CO2 concentration

The regulation and interaction of carbon and nitrogen metabolism is very important in plant metabolism [38]. Under higher CO2 concentrations, C3PMM predicted a decreased flux of nitrate ( NO-3 ) assimilation, as for photorespiration. In contrast, the flux of ammonium ( NH+4 ) assimilation increases with the CO2 concentration (Figure 6). The ratios of fluxes in photorespiration and ( NO-3 ) assimilation pathway are 0.728, 0.519 and 0.422 respectively under 550, 800 and 1000 ppm CO2 treatments, and the corresponding ratio of flux in ( NH+4 ) assimilation are 1.1, 1.177 and 1.212 under these three elevated CO2 conditions. In addition, the flux ranges of nitrate assimilation are [0, 10.4906], [0, 7.6363], [0, 5.4475], and [0, 4.4314] under 380, 550, 800, and 1000 ppm CO2 respectively, which is consistent with the flux shift of reactions in photorespiration pathway. While the flux range of ammonium assimilation are [0, 191.728], [0, 201.547], [0, 209.127], and [0, 212.642] under the four CO2 levels. It also supports the decreased flux of nitrate assimilation and increased flux of ammonium assimilation with the increasing CO2 concentrations.

This result is consistent with the experimental observations. Several studies by Bloom lab [1,39-41] demonstrated that higher CO2 inhibits nitrate assimilation in wheat and Arabidopsis, which would lead to lower organic N production. They also found similar responses in barley [42] and tomato [43]. Bunce showed that elevated CO2 inhibited ( NO-3 ) assimilation and carbohydrate translocation in soybean receiving both ( NH+4 ) and ( NO-3 ) as N sources [44].

Discussion

Coordination of plant primary metabolism

With the validated constraint-based model of C3 primary metabolism, we theoretically predicted the required changes in different reaction fluxes when CO2 concentration increases. Results showed that the predicted changes in the fluxes are consistent with the observed responses at the transcription level in the involved enzymes under elevated CO2 (Figure 5 and Table 3). This indicated that there is coordination between gene expression and metabolic flux in the leaf primary metabolism, to gain the maximal potential CO2 uptake rate. This finding is consistent with the long recognized high level of coordination in the prokaryotic Escherichia coli system [6]. Actually, the fundamental metabolic tasks are highly similar across divergent species, including nutrient uptake, the coordination of central carbon metabolism, the generation of energy, the supply of amino acids and protein synthesis [45]. As these tasks are universal across species, the input and output of specific regulatory modules are often similar for different organisms. The high level of coordination between reactions in the plant primary metabolism is also reflected in the relative small variance in the ratio between activities of different enzymes involved in the plant primary metabolism [46].

Furthermore, reactions in the same pathway show similar responses when the concentration of CO2 changes (Figure 5; Tables S1 and S2), which also reflect the coordination of plant primary metabolism. Having enzymes in the same pathway have the same level of changes helps to ensure a high resource use efficiency, particularly for use of nitrogen. This is because theoretically, for enzymes involved in a linear pathway, the most efficient way to distribute nitrogen between enzymes is to ensure that each enzyme has equal control coefficient over flux [47]. The coordinated expression at the transcription level might represent one of the regulatory cascades used by plants to achieve such coordination. It is reported in microbial metabolism, some other regulatory mechanisms, such as allosteric regulation by metabolite binding and post-translational protein modifications, were also shown to play important role for the coordination of metabolism [45]. We will consider the more complicated regulation to address the coordination of plant metabolism in future work.

Objective function of the mesophyll cell metabolism under light

Identifying objective function for each particular cell type is a major challenge when the constraint- based modeling approach is used in plant biological research. Analysis from this study demonstrated that maximization of CO2 uptake rate is the most possible objective function for C3 mesophyll cell under light. If the objective function is maximization of CO2 uptake rate, in general, the predicted changes in the flux under elevated CO2 are most consistent with the observed changes at the gene expression level (Table 3). Furthermore, analysis results suggested that the modern C3 plants have evolved mechanisms to respond to environmental perturbations in a manner that maximizes CO2 uptake, which is also consistent with the observed changes in photosynthesis, in particular Rubisco activity, under elevated CO2 [48,49].

Interaction between nitrogen and carbon metabolism under elevated CO2

With C3PMM and the identified objective function, we studied the interaction between carbon and nitrogen metabolism in a C3 mesophyll cell. The C3PMM predicted an inhibition of ( NO-3 ), but not ( NH+4 ), assimilation under elevated CO2 (Figure 6), which is consistent with the experimental observation in Arabidopsis and wheat [1,39-41], barley [42], tomato [43], and soybean [44]. One of the physiological mechanisms for this effect is that under elevated CO2 concentration, a decrease in photorespiration is expected, leading to a reduction in the NADH levels available to power nitrate assimilation. In contrast, no such phenomena in C4 plants, because the first carboxylation reaction in the C4 carbon fixation pathway generates adequate amounts of malate and NADH in the cytoplasm of mesophyll cells [50]. Another mechanism would be that elevated CO2 inhibits nitrite ( NO-2 ) transport from the cytoplasm into the chloroplast, the site where the subsequent conversion into amino acids occurs [41].

Previous studies have demonstrated that pine forests [51] or wetlands [52], where ( NH+4 ) is the dominant nitrogen form, showed a relatively large increase (25%) in net primary productivity under CO2 enrichment, whereas ecosystems in which ( NO-3 ) is the dominant nitrogen source, such as grasslands [53] or wheat fields [54], at standard fertilizer levels show declines in net primary productivity under CO2 enrichment. These results suggested that the decrease in photosynthesis and growth of plants after long exposures to elevated CO2, which is known as CO2 acclimation [49], might be due to decreased ( NO-3 ) assimilation, and the resultant decline in plant organic N [1,55].

Conclusion

In summary, plant primary metabolism is highly coordinated and therefore studying the interaction and coordination between different pathways are critical for the future rationale design of plant primary metabolism. Here we developed and used a constraint-based model of C3 Plant Primary Metabolism (C3PMM) to study the coordination between different pathways. We presented evidences supporting coordination of gene expression changes and metabolic flux changes under perturbed CO2 concentrations. Furthermore, maximizing photosynthetic CO2 uptake is the most possible objective function for a mesophyll cell. Flux balance analysis using this C3PMM suggested that the ( NO-3 ) assimilation was inhibited under elevated CO2, which suggested a mechanistic linkage between elevated CO2 and lack of CO2 fertilization effects in some ecosystems where the dominant source of nitrogen is ( NH+4 ). In conclusion, all the coordination of C3 plant primary metabolism at both transcriptome level and metabolism level, and the coordination between nitrate assimilation and photorespiration under elevated CO2 concentration, are beneficial for maximal CO2 uptake rate [56-60].

Acknowledgement

We thank Pinghua Li for providing us the complete microarray dataset of Arabidopsis thaliana under 550 vs 380ppm CO2 concentration. This work was supported by the Shanghai Natural Science Funding 16ZR1449700, the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry 15Z102050028, “SMC-SJTU Chen Xing” Program for Excellent Young Scholars 13X100010027, and open project from Key Laboratory of Synthetic Biology and Computational Biology.

Authors' Contributions

ZW designed the project and analyzed the C3 primary metabolism model. LYL reconstructed the model. LL conducted the analysis of metabolic network topology. JL compared the flux change and gene expression change. ZW and LL wrote the manuscript.