Potential metabolic mechanisms for inhibited chloroplast nitrogen assimilation under high CO2.

Improving photosynthesis is considered a major and feasible option to dramatically increase crop yield potential. Increased atmospheric CO2 concentration often stimulates both photosynthesis and crop yield, but decreases protein content in the main C3 cereal crops. This decreased protein content in crops constrains the benefits of elevated CO2 on crop yield and affects their nutritional value for humans. To support studies of photosynthetic nitrogen assimilation and its complex interaction with photosynthetic carbon metabolism for crop improvement, we developed a dynamic systems model of plant primary metabolism, which includes the Calvin-Benson cycle, the photorespiration pathway, starch synthesis, glycolysis-gluconeogenesis, the tricarboxylic acid cycle, and chloroplastic nitrogen assimilation. This model successfully captures responses of net photosynthetic CO2 uptake rate (A), respiration rate, and nitrogen assimilation rate to different irradiance and CO2 levels. We then used this model to predict inhibition of nitrogen assimilation under elevated CO2. The potential mechanisms underlying inhibited nitrogen assimilation under elevated CO2 were further explored with this model. Simulations suggest that enhancing the supply of α-ketoglutarate is a potential strategy to maintain high rates of nitrogen assimilation under elevated CO2. This model can be used as a heuristic tool to support research on interactions between photosynthesis, respiration, and nitrogen assimilation. It also provides a basic framework to support the design and engineering of C3 plant primary metabolism for enhanced photosynthetic efficiency and nitrogen assimilation in the coming high-CO2 world.

Introduction

With the continued increase in global population, and the increased urban population with improved diet, there is a great need to increase crop yield per unit land area. It is projected that, by 2050, the global food production needs to be doubled to meet the projected demand (Ray et al., 2013). Improvement of photosynthetic efficiency is regarded as a major feasible option to increase crop yield potential (Zhu et al., 2010). For C3 crops, one of the main strategies is inhibiting photorespiration by increasing the CO2 concentration around Rubisco (Long et al., 2006, 2006; Zhu et al., 2010). By the end of this century, the atmospheric CO2 concentration may reach 800 ppm, i.e. doubling the current atmospheric CO2 level (http://www.ipcc.ch). Elevated CO2 can influence both photosynthesis and other aspects of plant primary metabolism (Drake et al., 1997; Leakey et al., 2009; Dusenge et al., 2019). A large amount of data from free air CO2 enrichment (FACE) experiments show that elevated CO2 increases both photosynthesis and biomass, but the increase is less than the expectation (Ainsworth and Long, 2005; Long et al., 2006; Leakey et al., 2009; Kant et al., 2012). This is partly attributed to a decreased stomata conductance of C3 plants under elevated CO2 (Drake et al., 1997; Medlyn et al., 2001; Bernacchi et al., 2007), while limited nitrogen availability is another factor which constrains the increase in photosynthetic productivity (Reich et al., 2006; Terrer et al., 2019). Insufficient leaf nitrogen content can limit the sink development, such as grains, fruits, tubers, etc. (Ruiz-Vera et al., 2017), which causes feedback inhibition of photosynthetic gene expression (Paul and Foyer, 2001; Ainsworth and Long, 2005). Insufficient nitrogen also reduces the nitrogen investment into the photosynthetic apparatus, as both the Rubisco and chlorophyll content decrease under elevated CO2 concentration (Reich et al., 2006; Leakey et al., 2009; Kant et al., 2012; Dusenge et al., 2019). Decreased nitrogen assimilation and protein level in crops under elevated CO2 therefore may also become a serious threat to human nutrition (Myers et al., 2014).

Since the interactions of carbon metabolism and nitrogen metabolism are critical for the response of plants to elevated CO2 and the decreased nitrogen content in the plant may constrain the stimulation of photosynthetic productivity by elevated CO2 (Reich et al., 2006; Terrer et al., 2019), the impacts of elevated CO2 on nitrogen metabolism have been studied for many years. Long-term CO2 enrichment experiments (days to years) show that the decreased plant nitrogen content under elevated CO2 may be caused by the dilution effect, i.e. although the total nitrogen assimilation by the plant under elevated CO2 is increased, it is less than the enhancement of biomass, leading to lower nitrogen content on a mass basis (BassiriRad et al., 1997). In addition to this, elevated CO2 can decrease nitrogen absorption by the root system, nitrogen availability in the soil (Terrer et al., 2019), and nitrogen assimilation in the shoot (Smart et al., 1998; Bloom et al., 2010, 2012).

In illuminated C3 leaves, nitrate (NO3−) is converted to nitrite (NO2−) in the cytosol by nitrate reductase, with consumption of NADH. The NO2− is translocated into chloroplasts and further reduced to ammonium () catalyzed by nitrite reductase with consumption of FDH2. The in chloroplasts, either from reduction of NO3− or xylem absorption, is then fixed by a glutamine synthetase–glutamate synthetase (GS–GOGAT) complex with consumption of ATP, FDH2, and α-ketoglutarate (2-OG). The assimilation of NO3− and synthesis of are used to produce amino acids dependent on the energy, reductant, and carbon skeleton provided from carbon metabolism (Stitt et al., 2002; Nunes-Nesi et al., 2010; Foyer et al., 2011). The energy and reductant needed for carbon and nitrogen assimilation in chloroplasts are supported by the chloroplastic electron transport chain. Excess reductants and energy in the chloroplasts are exported to the cytosol by metabolite shuttles, e.g. the triose-photosphate shuttle (TPT) and the malate-oxoloacetate shuttle (OMT), or by photorespiration (Heineke et al., 1991; Igamberdiev et al., 2001; Stitt et al., 2002; Scheibe, 2004; Nunes-Nesi et al., 2010; Foyer et al., 2011).

For wheat (Triticum aestivum) and Arabidopsis thaliana exposed to elevated CO2 for multiple days, the chloroplastic nitrate assimilation rate is significantly inhibited (Bloom et al., 2010). In protoplasts treated with high CO2 conditions for minutes, the cytosolic reductant concentration decreases (Igamberdiev et al., 2001) because of the inhibited photorespiratory export of reducing power from chloroplasts to cytosol (Igamberdiev et al., 2001). Therefore, elevated CO2 is hypothesized to inhibit nitrogen assimilation in C3 leaves by decreasing the reductant required for reducing NO3− in the cytosol (Bloom et al., 2010, 2015). As a result, plants under elevated CO2 might preferentially absorb rather than NO3− (Bloom et al., 2010, 2012). This may partially explain the greater stimulation of Arabidopsis by elevated CO2 under treatment than that under NO3− treatment, as observed in growth chambers (Bloom et al., 2012). Different from this response, some C3 plants under elevated CO2 show no growth advantage under compared to NO3− (Andrews et al., 2019), while for others elevated CO2 increases NO3− absorption but decreases absorption (BassiriRad et al., 1997). On a whole plant basis, elevated CO2 can increase nitrogen assimilation (Andrews et al., 2019). To date, the mechanisms underlying the effects of elevated CO2 on chloroplastic nitrogen assimilation in C3 plants are still unresolved (Andrews et al., 2020; Bloom et al., 2020). These different responses might be related to the diversity of plant primary metabolism in C3 species, reflected in different enzymatic activities and metabolite levels (Arrivault et al., 2019).

Nitrogen assimlation is also closely linked to respiration, which makes the interactions between photosynthetic carbon and nitrogen assimilation even more complex. Mitochondrial respiration is inhibited in illuminated leaves (Tcherkez et al., 2008) and the tricarboxylic acid (TCA) cycle may even change from a closed mode to an open mode (Tcherkez et al., 2009; Sweetlove et al., 2010). This change might be partially attributed to the accumulation of certain metabolites under irradiance, e.g., the increased ratio of NADH/NAD+ in particular. Since NAD+ is the substrate of pyruvate dehydrogenase, isocitrate dehydrogenase, 2-OG dehydrogenase, and succinate dehydrogenase, there is competition for NAD+ and these enzymes are inhibited by NADH (Nunes-Nesi et al., 2013). The mode shift in the TCA cycle under illumination helps provide the carbon skeleton, in particular 2-OG, for nitrogen assimilation (Nunes-Nesi et al., 2010; Sweetlove et al., 2010). The 2-OG can not only be generated through the TCA cycle, but can also be generated by reactions catalyzed by aminotransferases (Chen and Gadal, 1990). For example, glutamate pyruvate transaminase (GPT) catalyzes the reversible conversion of glutamate (GLU) and pyruvate to alanine (ALA) and 2-OG; glutamate aspartate transaminase (GAT) catalyzes the reversible conversion of GLU and oxaloacetate (OAA) to aspartate (ASP) and 2-OG; and glutamate glyoxylate aminotransferase catalyzes the reversible conversion of GLU and glyoxylate to glycine (GLY) and 2-OG (Novitskaya et al., 2002; Sweetlove et al., 2010). Such diversity of plant primary metabolisms and the close interactions between photosynthesis, nitrogen assimilation, and respiration call for new approaches to study plant primary metabolism and to identify new options to manipulate plant primary metabolism for desired properties.

Dynamic systems models have been used to study complex metabolism networks and to identify engineering targets for crop improvements (Stitt et al., 2010; Baghalian et al., 2014; Zhu et al., 2016; Zhao et al., 2021). The major advantage of using systems models to guide crop engineering for increased efficiency is that this method can greatly decrease the time and labor cost associated with pre-screening the vast number of potential engineering options, as in the case of photosynthesis engineering (Long et al., 2015; Zhu et al., 2016, 2016). This method has been used successfully to show that increased recovery from photoprotective states can help improve photosynthetic efficiency and hence biomass production (Zhu et al., 2004; Kromdijk et al., 2016). So far, steady-state and dynamic systems models for both C3 and C4 photosynthesis have been developed, along with photosynthesis models at leaf scale and canopy scales (Farquhar et al., 1980; Pettersson and Rydepettersson, 1988; Laisk et al., 1989; Gross et al., 1991; Woodrow and Mott, 1993; Kirschbaum et al., 1997; Laisk and Edwards, 2000; Poolman et al., 2000; Laisk et al., 2006; Zhu et al., 2007, 2013; Song et al., 2013; Wang et al., 2014, 2020; Morales et al., 2018; Bellasio, 2019; Chang et al., 2019; Wu et al., 2019). For the photosynthetic metabolism models, the Farquhar, von Caemmerer and Berry (FvCB) model (1980) has been used widely to predict photosynthesis at the leaf level and beyond. This steady-state model was recently extended by integrating nitrogen assimilation with consideration of the interactions between photosynthesis, photorespiration and nitrogen assimilation (Busch et al., 2018). This extended model also predicted that increasing photorespiratory flux to supply carbon skeleton for nitrogen assimilation can increase net photosynthesis (Busch et al., 2018). However, the FvCB model and its improved versions (Farquhar et al., 1980; Harley and Sharkey, 1991; Busch et al., 2018) are stoichiometric models, which lack detailed description of enzyme kinetics and regulation of enzyme reactions by intermediate concentrations. As a result, these steady-state models cannot be used to predict the required genetic manipulations to improve photosynthetic efficiency and maintain nitrogen uptake.

Developing highly mechanistic kinetic models of photosynthesis and related reactions in the plant primary metabolism is desired to enable mechanistic prediction of the performance of plant primary metabolism under various genetic or environmental conditions (Fernie and Morgan, 2013; Zhu et al., 2016). Historically, dynamic systems models of photosynthesis with different mechanistic details were developed (Pettersson and Ryde-Pettersson, 1987, 1988; Poolman et al., 2000; Laisk et al., 2006; Zhu et al., 2007; Xin et al., 2015; Bellasio, 2019). These models have been used to study the mechanistic basis of different photosynthesis properties and identify new options to improve photosynthesis. However, in the models developed so far, the close interactions between photosynthesis, photorespiration, gluconeogenesis–glycolysis, the TCA cycle, and nitrogen assimilation (Stitt et al., 1984; Plaxton, 1996; Heldt, 2002; Riebeseel et al., 2010; Foyer et al., 2011; Timm et al., 2016; Heyneke and Fernie, 2018) are not considered.

In this article, we report the development of a dynamic systems model of plant primary metabolism, which can simultaneously simulate the fluxes of photosynthesis, glycolysis, the TCA cycle, and photo-nitrogen assimilation. This model realistically predicts the commonly observed responses of photosynthetic CO2 uptake rate (A), photorespiration (PR), day respiration (Rd), and nitrogen assimilation (NA) to variations of light and CO2 levels. Considering the interactions between photosynthetic carbon and nitrogen metabolism, we used this model to examine potential mechanisms underlying the observed responses in NA under elevated CO2. We envisage that this model or its extension will form a basic research tool to study interactions between photosynthesis, respiration and nitrogen assimilation. This model can also support identification of new options for engineering plant primary metabolism for increased photosynthetic efficiency, nitrogen assimilation or production of high-value products.

Results

Model parameterization

The abbreviations of metabolites and reactions were presented in Supplemental Tables S1 and S2, respectively. The model (Figure 1) was initially parameterized with data collected from literature (Supplemental Tables S3 and S4). The model was then used to simulate the net photosynthesis (A) under different CO2 and light levels (Supplemental Figure S1). The simulation results show that the model can realistically predict typical CO2 and light responses. Subsequently, Vmax was adjusted to minimize the difference between the modeled and measured response of A to CO2 (Supplemental Figure S4). The optimized Vmax of enzymes and transporters are listed in Supplemental Table S3.

Figure 1

The structure of the systems model of primary metabolism in a C3 mesophyll cell. In the C3 primary metabolism model, the cell was divided into four compartments including the chloroplast, the mitochondria, the peroxisome and cytosol. All abbreviations used in the figure were defined in the Supplemental Table S1.

Prediction of photosynthesis response to CO2 concentrations under different light and O2 conditions

We evaluated the capacity of the parameterized model to predict A under different light and O2 levels (Figure 2). The model realistically predicts the typical response curves of A to light (the AQ curve) and Ci (the A/Ci curve; experimental data in Figure 2 were the unpublished data gained from Dr Florian Busch in previous studies; Busch et al., 2018). Under saturating irradiance (photosynthetically active radiation [PPFD] = 1,800 μmol m−2 s−1), the modeled A shows a decrease as Ci rises above the saturating CO2 concentration, reflecting triose phosphate limitation at this CO2 range (Figure 2A). We then simulated the A/Ci curve under different O2 concentrations. Results show that, at a low Ci (<600 ppm), A under 2% O2 is always higher than that under 21% (Figure 2, A and D). This difference is due to the higher photorespiratory flux under a higher O2 concentration (Busch et al., 2018). However, under high Ci, A under 21% oxygen is consistently higher than that under 2.1% oxygen (Figure 2, A and D), reflecting triose phosphate limitation (Supplemental Fig. S3) under these conditions (Harley and Sharkey, 1991). We also used the model to simulate the A/Ci curves under different irradiances and AQ curves under different CO2 levels. The model successfully captured the behaviors of A under a couple of CO2 conditions and irradiances (Figure 2, B, C, E, and F).

Figure 2

Responses of net photosynthetic CO2 uptake (A) to environmental changes. (A and D) and (B and E) represent responses of A to CO2 concentration under different O2 and light conditions, respectively. (C and F) represent responses of A to photosynthetic photon flux density (PPFD) under different CO2 conditions. Solid lines (A, B, and C) represent simulation results; and dots and dotted lines (D, E, and F) represent experimental data, which were the unpublished data gained from Dr Florian Busch in previous studies (Busch et al., 2018). Data in parts (D, E, and F) are mean ± sd (n = 7 individual plants, the extra measurements at 400 ppm are double measurements on the same plants).

Prediction of the concentrations of metabolites of primary metabolism under different CO2 concentrations

Our model predicted the dynamic changes in concentrations of 96 metabolites with their subcellular locations fully considered (Figure 1). Among these metabolites, 25 of them have been measured previously in Arabidopsis (Arrivault et al., 2009). To enable a direct comparison between our model predicted and experimentally measured metabolite concentrations, we normalized our predicted metabolite concentrations to metabolite content on a leaf volume basis (Eq. 28). Figure 3 shows the comparison between measured and model-predicted metabolite concentrations under three different CO2 concentrations. Here, we describe the comparison between model simulation and experimental values for the metabolites involved in the Calvin–Benson cycle (CBC), starch/sucrose synthesis, and the glycolysis/TCA cycle.

Figure 3

The changes of modeled and measured metabolite profiles (Arrivault et al., 2009) under different CO2 conditions. The height of the bars represents the relative levels of metabolite in reference to that at the lowest CO2 condition. The blue, red, and dark yellow bars represent relative metabolite levels under CO2 levels of about 100, 200, and 500 ppm, respectively. The light density was 1,500 μmol m−2 s−1 in simulation and was 150 μmol m−2 s−1 in experimental treatment (A. thaliana). Simu, simulated metabolite levels; Meas, measured metabolite levels from Arrivault et al. (2009). All abbreviations used in the figure were defined in the Supplemental Table S1.

The model predicts increased concentrations for most of the metabolites in the CBC under elevated CO2 concentrations, which is generally consistent with the experimentally measured trends (Arrivault et al., 2009). Arrivault et al. (2009) showed that RuBP content was lower at the CO2 concentration of 200 ppm compared to [RuBP] under either 100 or 500 ppm (Figure 3). Our simulations show that [RuBP] initially increased then decreased with the increase of CO2 concentrations (Supplemental Figure S2). Furthermore, this pattern of RuBP concentration changes with Ci is maintained under different O2 and light conditions (Supplemental Figure S2). This simulated relation between [RuBP] and Ci is consistent with experimental data (Badger et al., 1984). The concentrations of PGA, SBP, and R5P in Arabidopsis were shown, by Arrivault et al. (2009), to be insensitive to CO2 changes. In our simulations, concentrations of all these substrates rose with increasing CO2 (Figure 3). Other studies carried out in C3 species showed that the PGA concentration gradually increases with CO2 (Badger et al., 1984; Leegood and von Caemmerer, 1994).

Most of the metabolites involved in starch and sucrose synthesis show increased concentrations with the increase in CO2 (Badger et al., 1984; Arrivault et al., 2009). This feature is generally captured by our model (Figure 3). One exception is the concentration of UDP glucose (UDPG). Our model predicts that [UDPG] gradually decreases with increasing CO2 concentration (Figure 3). However, in experiments, the UDPG concentration shows a gradual increase with increasing CO2 (Arrivault et al., 2009). The decrease in simulated [UDPG] results from a decrease in the model-simulated cytosolic ATP concentration (Figure 3).

PEP is an intermediate in glycolysis. The PEP concentration is known to increase with increasing CO2 levels (Leegood and von Caemmerer, 1994; Arrivault et al., 2009), which is successfully predicted by our model (Figure 3). In contrast, the concentration of 2-OG was predicted to decrease under elevated CO2 (Figure 3). Though the absolute values of these metabolites differ between measurements and simulations, the predicted patterns of changes in metabolite levels under different CO2 levels were largely consistent with experimental data (Arrivault et al., 2009).

The model also predicts the concentrations of the energy carriers in a cell, including NADP, NAD, and ATP. Consistent with measurements, the model predicts decreased NADP concentrations under high CO2 concentrations. Importantly, in the current model, there are also a number of metabolites for which the model prediction and experimental measurements differ. These metabolites include pyruvate, malate (MAL), succinic acid (SUCC), NAD, and ATP (Figure 3).

Predicted fluxes of photorespiration, respiration, and nitrogen assimilation under different light and CO2 concentrations

The fluxes of photorespiration under different CO2 and light levels were simulated using our model (Figure 4, A and D). Under saturated light (PPFD = 1500 μmol m−2 s−1) and ambient CO2 and O2 concentrations ([CO2] as 400 ppm, [O2] as 21%), the PR is 4.77 μmol m−2 s−1. Under the same conditions and model parameterization, the predicted A is 28.14 μmol m−2 s−1, resulting in a ratio of PR to A of 0.17, which is within the range of the ratio reported for C3 leaves (Perry et al., 1983). As expected, with an increase in Ci, the PR decreases (Figure 4A). Similarly, with an increase in PPFD, the model predicts an increase in PR (Figure 4D).

Figure 4

The predicted light and CO2 response curves of Rd, PR,NA. A–C represent CO2 response curves of PR, Rd, and NA, respectively. D–F represent light response curve of PR, Rd, and NA, respectively. The PPFD was 1,500 μmol m−2 s−1 in (A–C). The ambient CO2 was 400 ppm in (D–F). The dashed line in (A–C) corresponds to the ambient CO2 level of 400 ppm.

Our simulation shows that dark respiration (Rd) decreases with increasing irradiance. The decrease of Rd is rapid below 200 μmol m−2 s−1 and is slower at higher irradiances (Figure 4E). Under a saturating PPFD (1,500 μmol m−2 s−1) and ambient O2 concentration (21%), Rd is predicted as 0.78 μmol m−2 s−1. This value falls well within the range of reported Rd, i.e. 0.5–1.5 μmol m−2 s−1 for a C3 leaf (Tcherkez et al., 2017; Gauthier et al., 2018).

The effect of CO2 concentration on Rd remains ambiguous in the literature. Rd was reported to be enhanced (Leakey et al., 2009), repressed (Amthor et al., 1992; Tcherkez et al., 2008) or unaltered (Tcherkez et al., 2012) under elevated CO2 concentrations. Using our model, we predicted a decreased trend of Rd with increasing Ci (Figure 4B).

The predicted NA, including both and NO3− assimilation rates, is 1.28 μmol m−2 s−1 under a PPFD of 1,500 µmol m−2 s−1, a [CO2] of 400 ppm, and a [O2] of 21% (Table 1 and Figure 4C). This value is also within the range of the NA estimated previously (Busch et al., 2018). Tcherkez and Limami (2019) suggested that nitrogen assimilation rate is likely about 0.1 μmol m−2 s−1 (Tcherkez and Limami, 2019). Walker et al. (2014) estimated a nitrate assimilation rate in Arabidopsis according to previous studies (Rachmilevitch et al., 2004; Bloom et al., 2010) that is 0.07–0.6 μmol m−2 s−1 (Walker et al., 2014). At the present, there is no direct measurement of NA in leaves. The methods of measuring nitrogen assimilation flux in vivo need to be developed in the future. Under a saturating PPFD, NA is predicted to decrease until Ci reaches about 500–600 ppm, remaining constant at a higher Ci (Figure 4C). At an ambient CO2 level of 400 ppm, NA is predicted to increase rapidly under a low PPFD range. The rate then increases more slowly at higher irradiances (Figure 4F).

Table 1

Enhanced CO2 concentration inhibited nitrogen assimilation under different light levels. Nitrogen assimilation rate includes nitrate assimilation and ammonium assimilation. Simulations were done for two CO2 levels (ppm) and three light levels (μmol m−2 s−1). The percentage of changes of the assimilation rates of total nitrogen, nitrate, and ammonium under enhanced CO2 conditions compared with those under normal CO2 conditions were calculated. “−” means decrease. “+” means increase.

Irradiance	CO2	Nitrogen Assimilation		Nitrate Assimilation			Ammonium Assimilation
		Rate	Relative Change	Rate	Relative Change	Contribution	Rate	Relative Change	Contribution
μmol m-2 s-1	ppm	μmol m-2 s-1	%	μmol m-2 s-1	%	%	μmol m-2 s-1	%	%
300	400	1.17		0.44			0.73
300	800	1.10	−5.67	0.42	−5.14	34.02	0.69	−6.00	65.98
500	400	1.28		0.56			0.73
500	800	1.18	−8.19	0.54	−2.32	12.32	0.63	−12.69	87.68
1,500	400	1.28		0.81			0.47
1,500	800	1.16	−9.05	0.77	−4.27	29.84	0.39	−17.28	70.16

These nonlinear responses of NA to changes in either [CO2] or light levels imply a complex and mixed pattern of responses of nitrogen assimilation to different light and CO2 conditions. A number of earlier reports show that nitrogen assimilation is inhibited under elevated CO2 (Bloom et al., 2010; Dusenge et al., 2019). Under PPFDs of 300, 500, and 1,500 μmol m−2 s−1, the predicted rates of nitrate assimilation when pCO2 (ambient CO2 concentration) = 800 ppm were lower by 5.14%, 2.32%, and 4.27%, respectively, than rates simulated when pCO2 = 400 ppm. The predicted rates of ammonium assimilation under these three light levels when pCO2 = 800 ppm were 6.00%, 12.69%, and 17.28% lower compared to rates when pCO2 = 400 ppm. This comparison implies that the decline of nitrate assimilation is responsible for 34%, 12%, and 29% of the decrease of nitrogen assimilation under elevated CO2 under PPFDs of 300, 500, and 1,500 μmol m−2 s−1, respectively. On the other hand, the decline of ammonium assimilation is responsible for 66%, 88%, and 70% of the decreased nitrogen assimilation when pCO2 = 800 ppm, at irradiances of 300, 500, and 1,500 μmol m−2 s−1, respectively (Table 1).

The biochemical effects of elevated CO2 on chloroplastic nitrogen assimilation

After model validation, we then used the model to illustrate the biochemical mechanisms underlying the observed inhibition of nitrogen assimilation under elevated CO2. Simulations were conducted under the same conditions as those in Table 1. During the simulation, values of all enzymatic parameters were assumed constant and the nitrogen assimilation rate only depended on metabolite levels in this pathway. Under elevated CO2 (Figure 5) and different irradiances, the model predicted accumulated NO3−, which was observed earlier in C3 species (Smart et al., 1998; Bloom et al., 2010, 2012; Wujeska-Klause et al., 2019), a decrease in the cytosolic NADH/NAD+ ratio as also shown experimentally by (Igamberdiev et al., 2001), and a declined NO2− content in both cytosol and chloroplast (Figure 5). Interestingly, an increased concentration in the chloroplast, and the decreased ATP (Igamberdiev et al., 2001) and 2-OG (Novitskaya et al., 2002; Arrivault et al., 2009) contents under elevated CO2 were realistically predicted (Figure 5;Supplemental Figure S5).

Figure 5

The influence of elevated CO2 on nitrogen assimilation rate under different photosynthetic photon flux densities. The number of the x-axis represents PPFD during simulation with a unit of “μmol m−2 s−1”. The blue bars and red bars represent either concentrations of metabolites or particular ratio (e.g. NADH/NAD) values under CO2 concentrations of 400 ppm and 800 ppm, respectively. The “cyt” and “chl” in the subscripts represent their concentration (or ratio) in cytosol and chloroplast, respectively.

Why does the level of 2-OG decrease under elevated CO2? Simulations indicated that elevated CO2 increased the rates of 2-OG production through the TCA cycle, GPT, and glutamate OAA transaminase (Figure 6). However, the total increase in these rates was less than the decrease in the 2-OG generation rate as a result of the decreased photorespiration. Consequently, the level of 2-OG decreases with elevated CO2 (Figures 3 and 5;Supplemental Figure S5). Additionally, Novitskaya et al. also found that, in both wheat and potato (Solanum tuberosum), the ratio of GLY content to serine content (Gly/Ser) is decreased with enhanced CO2; the proportions of both ASP and alanine in total amino acids are increased with rising CO2 (Novitskaya et al., 2002). Our model predictions matched well with these experimental results (Supplemental Figure S7).

Figure 6

Model simulated metabolic fluxes of net 2-OG production to provide carbon skeletons as nitrogen assimilation in the C3 primary metabolism model, under three light levels and two CO2 conditions. The color in the bar graph represents the different reactions generating 2-OG, as shown in the left part of the figure. Substrates ICIT, glyoxylate (GOX), GLU, GLY, SER, PYR, ALA, ASP, OAA, and 2-OG represent iso-citrate, GOX, glutamate, GLY, serine, pyruvate, alanine, ASP, OAA, and 2-oxoglutarate, respectively. Enzymes GPT and GOT represent glutamate: pyruvate transaminase and glutamate: OAA transaminase, respectively.

Since redox levels play a crucial role in determining changes in nitrogen assimilation and the form of nitrogen assimilated, we further used the model to predict the level of reducing equivalents in chloroplasts under different light and CO2 levels. Under all simulated irradiances, the increased photosynthetic carboxylation due to elevated CO2 increased both the flux and the percentage of FDH2 consumed through CBC among all different FDH2 consumption pathways (Figure 7, A and B). Elevated CO2 caused chloroplastic FDH2 to rise slightly under low light while decreasing under high light (Figure 5). Model analysis shows that (Figure 7A), under elevated CO2 and sub-saturating irradiance (PPFD = 300 and 500 μmol m−2 s−1), the increased consumption rate of FDH2 by the CBC was less than the decrease in FDH2 consumption by OMT, NiR, and GS–GOGAT (Figure 7A). However, under a light level of 1500 μmol m−2 s−1, the increased consumption rate of FDH2 through CBC by elevated CO2 is greater than the total decrease of FDH2 consumptions by OMT, NiR, and GS–GOGAT. It implies that, in the chloroplast, competition for reducing equivalent between carbon assimilation and nitrogen assimilation is more intense under high PPFD. This conclusion should be treated with caution as the consumption of redox equivalents for the CBC activation was not considered in the current model.

Figure 7

Chloroplastic FDH2 consumptions through different metabolic pathways in under three light conditions and two CO2 levels. The left and right parts represent the fluxes and percentage of fluxes of FDH2 consumptions, respectively. The color of purple, yellow, orange, and blue represent the pathway of GS–GOGAT system, NO2-reduction (NiR), chloroplastic OAA/malate shuttle (OMT), and the CBC.

The predicted mass ratio of assimilated carbon to assimilated nitrogen in illuminated C3 leaves under different irradiances and CO2 conditions

Our model further simulates changes in A and NA under elevated CO2, which results in an increase in the mass ratio of carbon assimilation to nitrogen assimilation, here we name it as C/NMRA (Figure 8). This predicted C/NMRA ratio is in the range of 5–30, which is near the reported C/N ratio range of 8–90 in C3 plants (Kubiske and Pregitzer, 1996; Ruan et al., 2010). Additionally, the predicted C/NMRA was larger under high irradiances than that under low irradiance (Figure 8). These predictions are all consistent with previous observations that the nitrogen content in leaves is lower under higher irradiance or under elevated CO2 conditions (Kubiske and Pregitzer, 1996).

Figure 8

Changes in the ratio of fixed C to N (C/N) under different CO2 and light levels. The O2 concentration used in the simulation was 21%. The C/N was calculated by the equation C/N = (12*A)/(14*NA), where A and NA represent net photosynthetic CO2 uptake rate and nitrogen assimilation rate, respectively; 12 and 14 represent the molecular weights of C and N, respectively.

Discussion

A dynamic systems model of primary metabolism for C3 leaves

This article reports the development of a comprehensive dynamic systems model of C3 plant primary metabolism. This model is unique since it includes all the reactions of the CBC, photorespiration, starch and sucrose synthesis, glycolysis, the TCA cycle, nitrogen assimilation, and the steps involved in translocating metabolites between different cellular compartments (Figure 1). By accurately simulating the close interactions between photosynthetic carbon metabolism and nitrogen assimilation under different levels of irradiance and various CO2 concentrations, the model can be used to support a greater understanding of primary metabolism and the identification of potential targets for metabolic improvement with better parameterization.

It is important to acknowledge that many steady-state models of plant primary metabolism have been developed, which have contributed to insights on plant metabolism and its engineering (von Caemmerer, 2000; Farquhar et al., 2001; Poolman et al., 2004; Laisk et al., 2006; Sweetlove and Ratcliffe, 2011; Buckley et al., 2012; Sweetlove and Fernie, 2013; Cheung et al., 2014; Busch et al., 2018). One major advantage of these steady-state models is that they are easy to parameterize. By comparison, the dynamic systems model enables more flexible evaluation of the impacts of modifying environmental conditions, kinetic and regulatory properties of enzymes, etc., on a wide range of systems properties (Baghalian et al., 2014).

This model successfully reproduces the behaviors of C3 primary metabolism under different CO2 and irradiance levels

Any robust model of plant primary metabolism needs to be able to predict the responses of photosynthesis, photorespiration, day respiration, and nitrogen assimilation fluxes under different CO2 levels and irradiances. Our model simulates responses of A to intercellular CO2 concentrations (i.e. the A/C curve), and to PPFD (the A/Q curve), which are consistent with the experimental data (Figure 2). In particular, the model successfully predicts the reversed sensitivity of A to O2 under a saturating CO2 concentration (Figure 2, A and D) and shows that this is attributed to a phosphate (Pi) recycling limitation (Supplemental Figure S3; Sharkey, 1985, 1985; Busch et al., 2018). In addition to A, this model also captures features of the responses of Rd, PR, and NA to CO2 concentration (Figure 4, A–C) and irradiances (Figure 4, D–F). Under a PPFD of 1,500 μmol m−2 s−1, [CO2] of 400 ppm, and [O2] of 21%, the predicted values of PR (Figure 4A), Rd (Figure 4B), and NA (Figure 4C) are all within the ranges of experimentally measured or previously estimated values (Figure 4).

One important motivation to develop this model of plant primary metabolism is to improve predictions of metabolite concentrations involved in primary metabolism. This is important since, on one hand, metabolite concentra tions underlie the interactions between photosynthesis, respiration, and nitrogen assimilation; and on the other hand, a model capable of predicting dynamic changes of metabolites will result in a platform for designing new engineering options to increase photosynthetic efficiency, and the production of precursors for high-value products. The predicted (Figures 3 and 5;Supplemental Figures S2 and S5) responses of metabolite levels to CO2 have been compared with the experimental data (Badger et al., 1984; Leegood and von Caemmerer, 1994; Igamberdiev et al., 2001; Novitskaya et al., 2002; Arrivault et al., 2009). The responses of most (near 70%) metabolites to CO2 have been well predicted.

The current model failed to predict trends of concentration changes under higher CO2 concentrations for a number of metabolites, such as MAL, PYR, G1P, UDPG, SBP, and SUCC. The reason for this might be that some of these metabolites are distributed in multiple cellular compartments or used as substrates for multiple reactions, which are not yet included in the current model. The predicted contents of ATP and NAD in a leaf differ from the measurements (Figure 3), though the predicted ratios of NADH/NAD+ in the cytosol and of ATP/ADP in chloroplasts (Figure 5) agree with the measurements (Igamberdiev et al., 2001). Our model can accurately predict the fluxes of different metabolites (Figures 2 and 4). Here, we acknowledge that, for some metabolites, the model cannot predict their concentrations accurately, which warrants future efforts to improve the model. Measurement of metabolite concentrations in a compartment-specific manner is needed to gain a better understanding of the properties of primary metabolism and for better parameterization of the model.

Mechanisms underlying the responses of nitrogen assimilation to elevated CO2

To illustrate the use of our model for studying components of plant primary metabolism and their responses to climate change, we explored potential mechanisms underlying the response of NA to CO2 and irradiance. Since the values of parameters in our model are fixed, the predicted results represent the short-term (in minutes) effects of varying CO2 or irradiances on the C3 primary metabolism. We initially modeled the responses of NA to CO2 at saturating irradiance and the response of NA to irradiance at ambient CO2 and O2 levels. According to our simulations, NA decreased with CO2 until Ci was above 560 ppm (Figure 4C), and NA increased with irradiance until PPFD was ∼800 μmol m−2 s−1 (Figure 4F). Unfortunately, there were no direct measurements of nitrogen assimilation rate under different CO2 levels. Our predictions (Table 1) are consistent with a number of earlier observations where plants were treated under elevated CO2 for days to years (Gastal and Saugier, 1989; Bloom et al., 2010). This shows that the short-term response to elevated CO2 might, to some extent, reflects the metabolic properties of plants under long-term adaptation.

In the short-term, elevated CO2 decreases cytosolic NADH/NAD+ (Igamberdiev et al., 2001), which is consistent with the decreased nitrate reduction rate and hence a greater accumulation of NO3− in long-term elevated CO2 treatments (Smart et al., 1998; Bloom et al., 2010, 2012; Wujeska-Klause et al., 2019). Therefore, cytosolic redox power, rather than the availability of NO3−, limits nitrate assimilation under elevated CO2, as proposed earlier (Bloom et al., 2002, 2010; Bloom, 2015).

Similar to the NO3− content, an increased concentration in the chloroplasts was predicted under elevated CO2 (Figure 5). This difference suggests that elevated CO2 not only decreases nitrate reduction in the cytosol, but also inhibits ammonium fixation. The model with current parameterization predicted that, under elevated CO2, the inhibited ammonium fixation was responsible for >65% of the reduction of total nitrogen assimilation (Table 1). This inhibition, in short-term elevated CO2 conditions, can be attributed to the decreased concentrations of 2-OG (Novitskaya et al., 2002) and ATP in the chloroplast (Igamberdiev et al., 2001, 2001) under elevated CO2 (Figure 5;Supplemental Figure S5). Under three irradiances and two CO2 conditions, the modeled ATP concentrations in the chloroplast were all above 1 mM, while the predicted 2-OG concentrations in chloroplasts were about 6–20 μM (Figure 5). The Michaelis–Menten constant of GS for ATP is 0.6 mM (Acaster and Weitzman, 1985) and that of GOGAT for 2-OG is 0.07 mM (Lea et al., 1990). This comparison suggests that flux through the GS–GOGAT complex is more sensitive to decreasing [2-OG] than to decreasing [ATP].

Short-term elevated CO2 stimulates photosynthesis, resulting in production of more intermediates of carbon metabolism (Badger et al., 1984), potentially providing more carbon skeletons for nitrogen assimilation. However, both our model simulations (Figures 3, 5, 6;Supplemental Figure S5) and previous experimental observations (Novitskaya et al., 2002; Arrivault et al., 2009) show that the concentration of 2-OG decreases under elevated CO2 conditions. There is a positive relationship between photorespiration flux and 2-OG levels in leaves (Supplemental Figure S5; Novitskaya et al., 2002). Previous studies showed that GLY and serine can be formed through a transamination reaction with glutamate as a substrate in the photorespiratory pathway; a portion of the GLY and serine can be used to support synthesis of amino acids and proteins (Novitskaya et al., 2002). As a result, 2-OG can be generated from the photorespiratory pathway by transamination of glutamate (Scheible et al., 2000; Foyer and Noctor, 2002; Novitskaya et al., 2002).

As a substrate of the GS–GOGAT system, a decrease in 2-OG leads to a decreased rate of ammonia assimilation; indeed, maintenance of 2-OG levels to support nitrogen assimilation under elevated CO2 has previously been proposed (Stitt et al., 2002; Foyer et al., 2011). Interestingly, there are transcriptomics data for Arabidopsis, rice (Oryza sativa), and soybean (Glycine max) under long-term elevated CO2 treatment. In line with the discussion above, these data suggest that elevated CO2 increases the consumption of photosynthate by glycolysis and the TCA cycle to produce more energy and carbon skeletons (Li et al., 2008; Fukayama et al., 2009, 2011; Leakey et al., 2009). This result suggests that plants might have certain innate regulatory mechanisms to adjust their metabolism to maintain NA under long-term elevated CO2. More detailed molecular biology studies on it are needed to test this hypothesis and explore its implication for plant carbon and nitrogen status.

Our simulations also suggest that, with increasing irradiance, the preferred form of N for assimilation shifts from ammonium to nitrate (Table 1). This modeled change is a result of the increased reducing status of cytosol and chloroplast (Figure 5; Supplemental Figure S6), as indicated by the much higher NAD(P)H/NAD(P)+ in both compartments (Heineke et al., 1991; Igamberdiev and Gardeström, 2003). As discussed earlier, elevated CO2 inhibited both nitrate and ammonia assimilation (Table 1). The assimilation rate was more sensitive to elevated CO2 under all light levels (Table 1). This is because chloroplastic assimilation is more sensitive than cytosolic NO3− reduction to the decrease of redox equivalents. This difference in sensitivity is a result of the fact that the Michaelis constant of GOGAT for FDH2 (0.003 mM, Supplemental Table S3; Hirasawa et al., 1986) is much lower than that of chloroplastic nitrite reductase for FDH2 (0.006 mM, Supplemental Table S3; Mikami and Ida, 1989).

In this analysis, we mainly focus on the kinetic effects of changes in metabolite concentrations on the rates of primary metabolism. Specifically, we focus on the consequences of changes in metabolite levels on fluxes through NR, NiR, GOGAT, and FNR. Here, we emphasize that the response patterns of nitrogen assimilation to elevated CO2 in different species are reported earlier to be highly diverse. For example, contrary to our predictions, in some plants nitrate assimilation is more sensitive than ammonium assimilation to elevated CO2 (Bloom et al., 2002; Cheng et al., 2012; Rubio-Asensio and Bloom, 2017). But the declined absorption and increased NO3− absorption were also observed in C3 species (BassiriRad et al., 1997). Correspondingly, different species also show different changes in the protein content of biomass under elevated CO2. For example, in rice, where is the dominant form of assimilated nitrogen, the protein content of biomass decreases dramatically in a FACE experiment (Myers et al., 2014); in contrast, in many other species that use ammonia as the dominant form of assimilated nitrogen, much smaller changes in the protein content of biomass have been observed under elevated CO2 (Bloom et al., 2010). The mechanisms underlying these diverse responses are unknown, though it is likely that different kinetic properties of the enzymes related to nitrogen assimilation may play an important role. Moreover, our model simulations suggest that the reduction of nitrate assimilation and ammonium assimilation under elevated CO2 might be caused by decreased redox potential and decreased 2-OG content, respectively.

Critical research areas to help further improve the model

Finally, we emphasize that, in the current model simulations, the predicted responses of a few metabolites to CO2 and light levels differ from experimental measurements, indicating where further research is needed to improve this model and increase the accuracy of its predictions. First, we need a more accurate understanding of the metabolic structure of plant primary metabolism under different conditions and in different species. Though the basic structure of primary metabolism is well understood, this structure varies under different conditions. For example, the TCA cycle operates with different modes (Sweetlove et al., 2010; Michaeli et al., 2011; Batushansky et al., 2015); in particular, the oxidative pentose phosphate pathway and the glucose-6-phosphate shunt around the CBC can be activated in the dark and in low levels of irradiance (Sharkey and Weise, 2016; Preiser et al., 2019). The size and identity of different metabolically inactive pools in different cellular compartments also need to be clarified. Second, metabolite transport between cellular compartments needs to be better studied, including the localization, activity, kinetic properties, and rate laws governing various transporters in different organelle membranes. Third, recent data increasingly suggest large natural variations of metabolite concentrations in plants, reflecting diversity of primary metabolism in both C3 and C4 plants (Arrivault et al., 2019).

Furthermore, the current model was parameterized under saturated irradiance. However, many enzymes in plant primary metabolism, in particular those of the CBC, showed different activation statuses regulated by thioredoxin, pH, or Mg2+ (Woodrow and Berry, 1988). Multiple studies indicated that enzyme activities of the CBC were inactivated under dark and activated under illumination (Laing et al., 1981; Anderson et al., 1982; Wirtz et al., 1982; Kobza and Edwards, 1987; Hutchison et al., 2000) and the degree of their activations under illumination strongly depends on PPFD (Sassenrath-Cole and Pearcy, 1994; Sassenrath-Cole et al., 1994). To develop a more predictive model, the light regulation of enzyme activities should be integrated into this model (Borghi et al., 2019). Moreover, most of the metabolite regulations from previous publications and public databases were integrated into the current model. However, some intermediates, which may function as regulators for some enzymes of primary metabolism, have not been identified. This substantially weakens the predictive ability of a kinetic model (Hackett et al., 2016).

Lastly, the expression levels and the corresponding activities of different enzymes—including proteins involved in photosynthetic light reactions, enzymes, and metabolite transporters in plant primary metabolism—differ between different species, and between different conditions. These differences result from both changed gene expression levels (Weber and Fischer, 2007; Weber and Linka, 2011; Weber and Brautigam, 2013) and post-translational modifications (Wirtz et al., 1982; Sassenrath-Cole et al., 1994; Moorhead et al., 1999; Hodges et al., 2013; Shane et al., 2016). These factors can contribute to the diversity of metabolic fluxes, metabolite concentrations, and response behaviors in different environmental conditions and in different C3 species.

In summary, this study reports the development of a dynamic systems model of C3 primary metabolism that integrates photosynthetic carbon metabolism, respiration, and nitrogen metabolism. This model effectively summarizes huge amounts of physiological and biochemical knowledge accumulated during recent decades on C3 primary metabolism. Simulations with this model realistically capture many experimentally observed responses of C3 primary metabolism to different environmental conditions. Model simulations show that maintaining the 2-OG content in C3 leaves is critical for nitrogen assimilation under elevated CO2. This model provides a basic research platform to study C3 primary metabolism in silico. With better model parameterization, this model can also be used to identify potential engineering targets to manipulate for improved photosynthetic CO2 uptake and nitrogen assimilation.

Materials and methods

Model description

Development of the C3 primary metabolism model

The model of C3 primary metabolism in a mesophyll cell (Figure 1) was constructed by extending a previous model of photosynthetic carbon metabolism (Zhu et al., 2007) following a standard procedure for construction of dynamic systems models (Zhao et al., 2017). Here we describe the assumptions, rate equations, systems of ordinary differential equations, procedures used for model parameterization, and algorithms used to solve the systems of ordinary differential equations.

Assumptions of this model

The main assumptions used in the model include:

We only considered the interactions between photosynthetic carbon and nitrogen metabolism. The fluxes for lipid and nucleic acid synthesis in C3 leaves and their responses to changes in environmental conditions are not considered.

Since primary metabolism supplies the intermediates, ATP, and reducing equivalents for many cellular activities, almost all intermediates in primary metabolism can participate in multiple reactions in addition to their role in primary metabolism. In the current model, we assumed some metabolite sinks, i.e. malate, sucrose, starch, glutamate, ASP, alanine, GLY, serine, NADH, and ATP. We assumed Michaelis–Menten kinetics to describe the consumption or storage of these intermediates (Barenholz et al., 2017). For some metabolite sinks, the Hill constant (we empirically set it as 2 in some equations) was used since it produces a sigmoid curve. We chose this since such a response curve can dramatically decrease the sink consumption rate when the sink metabolite concentration is small, and dramatically increase the sink consumption rate when the sink metabolite concentration is large. In other words, such a sigmoid response pattern promotes a fast response of sink consumption rate to the change of sink metabolite concentration, which helps avoid large changes in metabolite concentrations and increases the robustness of such a complex metabolism model (Alon, 2006).

We assumed that all sink metabolites have no feedback regulation on their generations.

Following earlier models (Pettersson and Rydepettersson, 1988; Zhu et al., 2007, 2013; Wang et al., 2014), the current model also assumes that some limiting resources are maintained at constant levels. These limiting resources are the total concentrations of ATP, ADP, and ADPG in chloroplasts; FDH2 and FD in chloroplasts; NADP and NADPH in chloroplasts; ATP and ADP in cytosol; UTP, UDP, and UDPG in cytosol; NADPH and NADP in cytosol; NADH and NAD in cytosol; ATP and ADP in mitochondria; and NADH and NAD in mitochondria. The total concentration of Pi (including the inorganic Pi free and organic Pi in the metabolites considered in our model) in chloroplasts, cytosol, or mitochondria is also maintained at a constant level. The total concentration of Pi in each compartment is calculated as the sum of total Pi in different metabolites. The concentration of Pi in each metabolite is calculated as the product of the metabolite concentration and a Pi coefficient for this metabolite. If there is no Pi group in one metabolite, such as malate and most metabolites in the TCA cycle and nitrogen assimilation pathway, its Pi coefficient is 0; if there is only one Pi group in one metabolite, its Pi coefficient is 1, such as PGA and G6P and most intermediates in the CBC and glycolysis; if there are two Pi groups in one metabolite, its Pi coefficient is 2, such as RuBP and FBP. It should be noted that in our model and in previous models (Pettersson and Rydepettersson, 1988; Zhu et al., 2007, 2013; Wang et al., 2014), the Pi coefficient for ATP and UTP is 1 instead of 3 since only the conversions between ATP and ADP and between UTP and UDP were considered. During these conversions, only one Pi is involved in the reactions. For NADP and NADPH, their Pi coefficients are both 0 since we only considered the oxidation-reduction reaction between NADPH and NADP. So the Pi coefficient in our model does not represent the number of Pi groups in a metabolite, rather it represents the number of Pi groups that can be transferred between compounds in biochemical reactions in our model.

We defined the net CO2 assimilation rate in this model to be the difference between the sum of fluxes for all reactions consuming CO2 (reactions catalyzed by Rubisco and PEPC) and the sum of fluxes for all reactions releasing CO2 (reactions catalyzed by GDC complex, PDC complex, NAD-ME, and reactions of the TCA cycle).

The model did not consider storage of either nitrate or ammonium in the vacuole or NH3 gas emission.

The reactions catalyzed by the GS–GOGAT enzymatic complex were simplified as one reaction in this model. Such a simplification for this autocatalytic cycle dramatically increases the stability and robustness of this model (Barenholz et al., 2017). We assumed that there was no Gln exportation from primary metabolism.

The reactions of the H2O2–O2 cycle in the photorespiration pathway were not considered.

To enable the model to be used to guide rice metabolic engineering, we include the PEPC catalyzed PEP carboxylation in chloroplasts in the current model. This reaction has been reported in rice from independent studies (Masumoto et al., 2010; Muramatsu et al., 2015; Covshoff et al., 2016) where the decreased expression of chloroplastic PEPC resulted in stunted vegetative growth and dramatic changes in primary metabolism, especially under ammonium treatment (Masumoto et al., 2010). Its enzymatic activity in rice chloroplast and kinetic properties have been systematically studied (Masumoto et al., 2010; Muramatsu et al., 2015). Although the enzymatic activity of enolase and phosphoglycerate mutase (PGM) has not been studied in other C3 plants, the protein ENOLASE1 has been found in the chloroplast (www.uniprot.org/uniprot/Q9C9C4), and ENOLASE1 is predicted to have a protein–protein interaction with PGM (https://string-db.org/network/3702.AT3G12780.1). So, the reactions catalyzed by them to generate PEP from PGA in the chloroplast are constructed in this model. The fluxes of this pathway in C3 chloroplasts of other plants need to be confirmed by experiments in the future.

Rate equations

Metabolic reactions

The basic procedure used to develop rate equations for most of the reactions in the CBC, photorespiratory pathways, starch and sucrose synthesis, TCA cycle, glycolysis, and nitrogen assimilation follows Zhu et al. (2007). Here we first briefly describe the equations used to represent metabolic reactions. The full list of equations is given in Supplemental Appendix A1. The full list of abbreviations and the values used in the model are listed in Supplemental Appendix A2.

Most metabolic reactions included in the model can be described by standard Michaelis–Menten rate equations. These equations are grouped into reversible and irreversible reactions. The reversible reactions are described as:

Type 1: A ←→ B

Its rate equation is formulated as:

Type 2: A + B ←→ P + Q Its rate equation is formulated as:

In Equations 1 and 2, [A], [B], [P], and [Q] represent the concentrations of A, B, P and Q, respectively. Vmax represents the maximal velocity of an enzyme. K represents the equilibrium constant of this reaction. represent the Michaelis–Menten constants of the enzyme for metabolites A, B, P, and Q, respectively.

Rate equations for irreversible reactions are formulated as:

Enzymes can be either inhibited or activated by metabolites. An inhibitor can competitively inhibit an enzyme by modifying this inhibitor’s Michaelis–Menten constant (K) for the enzyme, leading to an increased apparent Michaelis–Menten constant (K):

This inhibitor can also noncompetitively inhibit an enzyme by modifying the V of this enzyme, leading to a decreased apparent Vmax ():

In Equations 5 and 6, the K represents the inhibition constant of metabolite X to the enzyme, and [X] represents the concentration of X.

Similarly, an activator can activate an enzyme by decreasing its apparent K ():

An activator can also directly increase the apparent Vmax for the enzyme :

In Equations 7 and 8, [X] represents the concentration of an activator. K represents the activation constant of metabolite X.

Following we describe reactions that cannot be described by Michaelis–Menten equations and are not included in Zhu et al. (2007).

Light reactions

The light reactions are described by the following equations (Eqs 9–15): represents the photon flux density of incident PPFD absorbed by photosystem II (PSII). I represents the incident PPFD. abs represents the proportion of PPFD absorbed by leaves; here it is assumed as 0.85 (Evans, 1987). f represents the proportion of light that cannot be used for photochemistry (assumed to be 0.15; Evans, 1987). 0.5 represents the proportion of light absorbed by PSII. J represents the photosynthetic electron transport rate. Jmax represents the maximal photosynthetic electron transport rate. represents the empirical curvature factor of the electron transport rate-irradiance curve (assumed to be 0.7; Evans, 1989).

Equations 11 and 12 represent equations used to calculate the rate of ATP synthesis by the photosynthetic electron transfer chain. represents the maximal rate of ATP synthase. represents the stoichiometric ratio of ATP production to electron transport (default value is 1 which is estimated from Caemmerer [2000]).

Equations 13 and 14 calculate the rate with which reduced ferrodoxin (FDH2) is produced through the photosynthetic electron transfer chain. represents the maximal rate of ferredoxin reductase in the chloroplast. is the stoichiometric number of Fd molecules reduced to FDH2 when 1 e- passes through the photosynthetic electron transfer chain (Carrillo and Ceccarelli, 2003). The default value of FDHyield is 0.5 which is estimated by referencing (Caemmerer, 2000).

Mitochondrial electron transport chain coupled with oxidation phosphorylation

The mitochondrial electron transport chain is simplified to be four reactions converting NADH (Reaction 1) and NADPH (Reaction 2) to NAD and NADP, respectively, in the cytosol, and converting NADH (Reaction 1) and FADH (Reaction 5) to NAD and FAD, respectively, in mitochondria. These oxidations of NADPH, NADH and FDH are coupled with oxidative phosphorylation to produce ATP. The oxidation of 1 molecule of mitochondrial NADH through complex I, III and IV transports 2e− to O2. It simultaneously exports 10 H+ from the mitochondrial matrix to the intermembrane space. The oxidations of cytosolic NAD(P)H by the NAD(P)H dehydrogenases or mitochondrial FADH by complex II each export 6 H+ from the mitochondria matrix to the intermembrane space while they simultaneously transport 2e− to O2 (O'Leary, 2016). This model does not include the alternative oxidase or other unique plant electron transport processes. This model assumes that 4H+ passing through ATP synthetase are needed to synthesize 1 ATP molecule according to the previous study (Wikström and Hummer, 2012; Hahn et al., 2018). Consequently, during oxidative phosphorylation in the mitochondria, oxidation of 1 NADH produces 2.5 molecules of ATP and oxidation of FADH produces 1.5 molecules. In the cytosol, NAD(P)H oxidation produces only 1.5 ATP molecules.

The reactions of NAD(P)H and FADH oxidation coupled with ATP synthesis (i.e. oxidation phosphorylation) are described as Michaelis–Menten equations as follows (Eqs 15–18):

Reaction 1:

Reaction 2:

Reaction 3:

Reaction 4:

In Equations 15–18, and represent the maximum velocities of NADH dehydrogenase and NADPH dehydrogenase, respectively, on the side corresponding to intermembrane space. represent the maximum velocities of mitochondrial complex I and complex II, respectively. and Pi represent the cytosolic concentrations of NADPH, NADH, ADP and Pi, respectively. represent the concentrations of mitochondrial FADH, NADH, ADP, and Pi, respectively; represents the Michaelis–Menten constant of cytosolic side NADH dehydrogenase for NADH. and represent the Michaelis–Menten constants of ATP synthase for ADP and Pi, respectively; and represent the Michaelis–Menten constants of complex I for NADH and of complex II for FADH, respectively.

Transporters

The metabolite transport process is described as a single-substrate reversible reaction with Michaelis–Menten equations (Eq. 19). So, for a metabolite M to be transported from outside (Mout) to inside (Min):

Min ←→Mout

V max represents the maximal velocity of the transporter, KM represents the Michaelis–Menten constant of the transported substrate M, and KE represents the equilibrium constant of this transporter, which is estimated based on the ratio of the “outside” concentration to the “inside” concentration from previous studies. Here, we describe the rate equations for different transporters involved in this study.

DiT1/2

The transporters DiT1 and DiT2 translocate 2-OG from cytosol to chloroplast and glutamate from chloroplast to cytosol, respectively. Their equations are adapted from Equation 19, as follows, to assure the robustness of model:

TPT

TPT couples export of DHAP, GAP or PGA from chloroplast to cytosol with an import of Pi. Their reactions were modified from the equation used by Zhu et al (2007):

ATP counter-exchange transporter on mitochondrial membrane

The ATP counter-exchange transporter on the mitochondrial inner membrane functions to export mitochondrial ATP to the cytosol and simultaneously to import ADP and Pi from the cytosol to the mitochondria. The equation describing this reaction has been developed in previous models (Bohnensack, 1982; Cortassa et al., 2003) and the constant 0.0203 was estimated from (Bohnensack, 1982; Cortassa et al., 2003):

Metabolite sinks

Some metabolites in the model can be consumed through either biochemical reactions or metabolite translocations. For example, sucrose, glutamate, and ASP can be imported into phloem and transported to other tissues to support their growth (Winter et al., 1992). Malate, alanine, GLY, and serine can be consumed by biochemical pathways or stored in the cell (Martinoia et al., 1985; Martinoia and Rentsch, 1994; Timm et al., 2013). In this model, the sinks of metabolite (including: malate, sucrose, starch, glutamate, ASP, alanine, GLY, serine, NADH, and ATP) are described with Michaelis–Menten equations (Zhu et al., 2007; Wang et al., 2014): hc is the Hill constant (Doran, 2013). To test the robustness of these parameters on predicted metabolic fluxes, we performed sensitivity analyses for them. The parameters (Vmax, K, and Hill coefficient) in the sink descriptions were perturbated through varying the default value by −10%, −5%, −3%, −1%, 1%, 3%, 5%, and 10% (Supplemental Figure S8).

Systems of ordinary differential equations and linear algebraic equations

The differential equations were derived following Zhu et al (2007). We similarly developed a set of algebraic equations describing mass balance relationships. The full list of mass balance equations is in Supplemental Appendix A3.

During our initial model parameterization, kinetic parameters, e.g. K, Ka, K, and K, were collected from the literature. Whenever possible, we used Vmax values obtained from sunlit leaves of C3 plants. The number of Vmax values in this model is more than the measured ones; these unmeasured Vmax(s) were first assumed to make the corresponding fluxes of their pathway in previously reported physiology ranges. The sources for each parameter used in the model are listed in Supplemental Appendix A1. The full list of enzyme activities is provided in Supplemental Appendix A2.

Since the protein contents of enzymes and enzyme activation status vary dependent on plant species, growth conditions, and developmental stages, we include an extra data-fitting step to improve model parameterization. In this step, we fitted the model to A/Ci curves measured under different nitrogen conditions (Busch et al., 2018). A genetic algorithm was used to optimize Vmax in the model by minimizing the differences between modeled and measured photosynthetic rates under different CO2 conditions (Eq. 28). The algorithm allows Vmax to vary in a fold range of 0.2–5. and represent the modeled value and the average value of measured A, respectively. σ represents the standard deviation of measured A.

Following Zhu et al. (2007), in this model, the initial default metabolite concentrations are based on reported metabolite concentrations for mature illuminated leaves (von Caemmerer and Edmondson, 1986; Winter et al., 1993, 1994; Leidreiter et al., 1995; Arrivault et al., 2009).

Numerical simulation experiments

The system of ordinary differential equations and the linear algebraic equations describing the model of plant primary metabolism are solved by using ode15s of the MATLAB package (MathWorks, ver 2012, Natick, MA, USA). Metabolite concentrations and rates of reactions or processes are the outputs of the model simulation runs.

Calculating the rate of A, PR, Rd, and NA

v RC represents the Rubisco carboxylation rate; vPEPC represents the CO2 fixed by PEPC; PR represents the rate of CO2 released by GDC (vGDC); Rd represents the total rate of CO2 released by PDH (vPDH), IDH (vIDH, including cytosolic and mitochondrial types), OGDH (vOGDH) and NAD-ME in mitochondria (vME). NA represents the sum of the rates of (vammonium_ass) and NO3− assimilation rates (vnitrite_ass). Only these five amino acids are synthesized in the current model since other amino acids can be synthesized from them as summarized (Heldt and Piechulla, 2011). They were chosen since they are directly linked to photosynthetic carbon metabolism and closely influence the C–N interaction.

Metabolite concentrations from the subcellular level to the leaf level

In this model, we predicted the metabolite concentrations in different subcellular compartments. To enable direct comparison of the predicted metabolite concentrations with measured whole-leaf metabolite concentrations (Arrivault et al., 2009), the predicted metabolite concentrations are converted to concentrations on a leaf volume basis by the following equation:

[M] represents the concentration of metabolite ‘M’ and the unit is mM. The subscripts chloroplast, cytosol, mitochondria, and leaf represent M concentration in chloroplast, cytosol, mitochondria, and leaf, respectively. , represent the area-based volume of chloroplast (30 mL·m−2), cytosol (30 mL·m−2), and mitochondria (4 mL·m−2), respectively.

Supplemental data

The following materials are available in the online version of this article.

The predicted net photosynthetic CO2 uptake rate (A) under different light and CO2 levels. This prediction was conducted before parameterization by a genetic algorithm.

Predicted concentrations of metabolites of the CBC under different conditions.

Predicted chloroplastic Pi (left) and ATP (right) concentrations under different intercellular CO2 levels (Ci) at either low or normal O2 concentrations.

The fold changes at different Vmax values in the model after parameterization using a genetic algorithm.

The predicted leaf 2-OG concentrations under different CO2 concentrations.

The responses of cytosolic (A and B) and chloroplastic (C) redox statuses under different PPFD (μmol m−2 s−1).

The predicted ratio of Gly/Ser (A), ASP concentration (B), and alanine (C) under different CO2 concentrations.

Sensitivity analyses for parameters used in the rate equation for sink capacity showed robust predictions of C3 primary metabolism fluxes.

Abbreviations have been used in the C3 primary metabolism model and in the manuscript.

The reactions or metabolic processes in the C3 primary metabolism model.

The corresponding maximum velocity of each enzyme, transporter of process in the C3 primary metabolism model.

The parameters and their values of Michaelis–Menten constants (KM), inhibition constants (KI), activation constants (KA), and equilibrium constants (KE) in the primary metabolism model.

Rate equations in the model.

Parameters and reactions in the model.

Differential equations in the model.

Click here for additional data file.