Towards a dynamic photosynthesis model to guide yield improvement in C4 crops

Summary The most productive C4 food and biofuel crops, such as Saccharum officinarum (sugarcane), Sorghum bicolor (sorghum) and Zea mays (maize), all use NADP‐ME‐type C4 photosynthesis. Despite high productivities, these crops fall well short of the theoretical maximum solar conversion efficiency of 6%. Understanding the basis of these inefficiencies is key for bioengineering and breeding strategies to increase the sustainable productivity of these major C4 crops. Photosynthesis is studied predominantly at steady state in saturating light. In field stands of these crops light is continually changing, and often with rapid fluctuations. Although light may change in a second, the adjustment of photosynthesis may take many minutes, leading to inefficiencies. We measured the rates of CO2 uptake and stomatal conductance of maize, sorghum and sugarcane under fluctuating light regimes. The gas exchange results were combined with a new dynamic photosynthesis model to infer the limiting factors under non‐steady‐state conditions. The dynamic photosynthesis model was developed from an existing C4 metabolic model for maize and extended to include: (i) post‐translational regulation of key photosynthetic enzymes and their temperature responses; (ii) dynamic stomatal conductance; and (iii) leaf energy balance. Testing the model outputs against measured rates of leaf CO2 uptake and stomatal conductance in the three C4 crops indicated that Rubisco activase, the pyruvate phosphate dikinase regulatory protein and stomatal conductance are the major limitations to the efficiency of NADP‐ME‐type C4 photosynthesis during dark‐to‐high light transitions. We propose that the level of influence of these limiting factors make them targets for bioengineering the improved photosynthetic efficiency of these key crops.


INTRODUCTION
Increasing crop genetic yield potential during the Green Revolution proved critical both to global food security and reducing the land-use change that would otherwise have been necessary to support the growing world population. The biggest contribution from this period was from new genotypes of major grain crops, which had greatly improved harvest indices and, coupled with improved agronomy, provided increased yield UNI-CEF, 2018, 2020). However, after rapid increases in yield over the 45 years from 1960, the growth in productivity of the world's major crops is stagnating (Ray et al., 2012). With a forecast 60% increase in global demand for primary foodstuffs by mid-century there is an urgent need to find new means to sustainably increase the productivity of these key crops (FAO, IFAD and UNICEF, 2020).
In parallel, a second global challenge is how to provide sustainable sources of energy and bioproducts to meet increasing societal pressures to achieve zero net glasshouse gas emissions. Displacing 55-70% of the petroleum demand in the USA with cellulosic biofuels would require about 10 9 Mg of biomass (Langholtz et al., 2016;Robertson et al., 2017). Saccharum officinarum (sugarcane) and Zea mays (maize) are the largest current sources of biofuels, whereas Sorghum bicolor (sweet and fiber sorghum) and the sugarcane-derivative fiber crop, energycane, have emerged as major potential bioenergy and bioproduct feedstocks (Crow et al., 2020;Gautam et al., 2020;Jaiswal et al., 2017;Long et al., 2015;Parajuli et al., 2020). All belong to the monophyletic grass tribe Andropogoneae and all predominantly use the C4 NADP-ME metabolic pathway (Bianconi et al., 2019). Increasing crop photosynthetic efficiency is one means to meet the need to increase the yields of food, biofuels and bioproduct crops and so avoid a need to bring more land into agriculture.
The yield potential of a given genotype at a given location is the product of the incident photosynthetically active radiation over the growing season, the efficiency of the crop in intercepting that radiation (ϵ i ), the efficiency of the conversion of intercepted radiation into plant mass (ϵ c ) and the efficiency of partitioning that mass into the harvested product (ϵ p ), also termed harvest index (Zhu et al., 2010); see Table 1 for a list of abbreviations, their definitions and units. Plant breeding has optimized ϵ i and ϵ p to the point where there is little opportunity for further improvement of these efficiencies in the major crops (Long et al., 2019;Zhu et al., 2010). Harvest indices are now as high as 0.6 in modern crop cultivars, with little scope for further improvement in this trait (Evans, 1997;Murchie et al., 2009;Zhu et al., 2010). Under optimal conditions the ϵ p of maize hybrids is about 0.52, which has remained unchanged for the last two decades (Di Matteo et al., 2016), whereas sugarcane has an already high ϵ p of 0.8, as the majority of the plant is the harvested product (i.e. stem) (Waclawovsky et al., 2010). By contrast, the third factor ϵ c (also known as lightuse efficiency) of C3 and C4 crops, governed by photosynthesis, falls well below its theoretical maximum (Zhu et al., 2008). Theoretical analysis and genetic engineering have shown considerable potential to increase photosynthetic efficiency in both C3 and C4 crops (Kromdijk et al., 2016;Long et al., 2019;López-Calcagno et al., 2020;Murchie et al., 2009;Salesse-Smith et al., 2018;South et al., 2019;Yoon et al., 2020). Although increasing the stress tolerance of crops is another important route to increasing productivity, experience shows that increasing genetic yield potential can address both of these factors, thereby increasing the average yields achieved under both optimal and stress conditions (Wu et al., 2019). For example, a detailed analysis of progressive gains in yield potential through soybean breeding have resulted in achieved yield increases in years with both good and suboptimal production conditions (Koester et al., 2014).
Mathematical modeling of the photosynthetic process and the application of optimization methods have proven valuable in identifying targets for the bioengineering of increased efficiency and sustainability (Wang et al., 2014;Zhu et al., 2004Zhu et al., , 2007. They have led to proof-of-concept improvements in productivity in replicated field trials with genetically engineered Nicotiana tabacum (tobacco) (Kromdijk et al., 2016;López-Calcagno et al., 2020;South et al., 2019). However, these models have focused largely on the steady state and often on light-saturated conditions, as have most measurements and analyses of the limitations of leaf photosynthesis.
In the field, however, leaves are rarely in steady state, and are subject to frequent fluctuations in light intensity. There has been a growing awareness of the need to address photosynthetic efficiency in fluctuating light conditions, as well as under constant light (Acevedo-Siaca et al., 2020b;De Souza et al., 2020;Deans et al., 2019;Hubbart et al., 2012;McAusland and Murchie, 2020;McAusland et al., 2016;Murchie and Ruban, 2020).
A great deal of progress has been made towards understanding the dynamic response to light in C3 plants. The non-steady-state photosynthetic rate of C3 plants under fluctuating light is mainly affected by the speed of the following responses: the relaxation of non-photochemical quenching; the activation and de-activation of Rubisco; the activation of ribulose 1,5-bisphosphate (RuBP) regeneration enzymes; metabolite pool sizes; mesophyll conductance and stomatal conductance (Deans et al., 2019;Kaiser et al., 2016;Mott and Woodrow, 2000;Pearcy, 1994;Slattery et al., 2018). However, the major limitations vary with species (Acevedo-Siaca et al., 2020;De Souza et al., 2020;McAusland et al., 2016;Taylor and Long, 2017).
Dynamic photosynthetic models infer that faster responses of photosynthesis to light fluctuations within canopies could increase the productivity of C3 crops by 13-32% (Taylor and Long, 2017;Wang et al., 2020;Zhu et al., 2004). A bioengineered increase in the speed of the relaxation of non-photochemical quenching, upon transition from high light to low light, significantly increased the quantum yield of photosynthesis in tobacco grown under fluctuating light, and also resulted in an approximately 20% increase in productivity in replicated field trials (Kromdijk et al., 2016).
Few experimental studies of photosynthesis under fluctuating light conditions in C4 plants have been published, however. The loss of productivity in fluctuating light, relative to constant light, in two C4 species (Amaranthus caudatus and Setaria macrostachya) was greater than in two C3 species (Celosia argentea and Triticum aestivum), implying a greater impact on C4 plants (Kubásek et al., 2013). Aside from this experiment, two early sets of gas-exchange measurements suggested that stomatal conductance is not a limiting factor during photosynthetic induction in maize (Furbank and Walker, 1985;Usuda and Edwards, 1984).
Energy-use efficiency limitations of C4 photosynthesis under steady-state conditions have been analyzed via a number of empirical analyses and biochemical models (Bellasio and Griffiths, 2014;Laisk and Edwards, 2000;Wang et al., 2014a,b;Yin and Struik, 2018;Yin and Struik, 2021). However, no mechanistic modeling studies have analyzed the limitations of non-steady-state C4 photosynthesis. Here, we developed a dynamic model to predict the potential limitations within C4 photosynthesis in fluctuating light to suggest feasible targets to improve the energyuse efficiency of C4 crops. As the major food and fiber C4 crops maize, sorghum, sugarcane and Miscanthus predominantly use the NADP-ME form of C4 photosynthesis, we previously developed a kinetic metabolic model of this form of C4 photosynthesis. It represented all discrete metabolic reactions of photosynthesis within the different cellular compartments of a C4 leaf, and transfer between the compartments, as a system of ordinary differential equations (ODEs). Each individual reaction and transfer between compartments is described using either enzyme kinetics or metabolite diffusion kinetics (Wang et al., 2014). This metabolic model was used to predict the limitations of C4 photosynthesis under steady-state lighting at 25°C.
To predict limitations under non-steady-state conditions, this metabolic model was extended to capture the key factors affecting non-steady-state photosynthesis during transitions from low light to high light, and vice versa. Specifically, we included the post-translational regulation of key photosynthetic enzymes, temperature responses of the enzyme activities, dynamic stomatal conductance and leaf energy balance. We then parameterized the model using gas-exchange data for the three most widely grown NADP-ME C4 crops: maize, sorghum and sugarcane. From this, the major factors limiting photosynthesis during transitions from the dark to high light for these crops were predicted to be Rubisco activase (Rca), the pyruvate, phosphate dikinase (PPDK) regulatory protein and stomatal conductance.

Factors influencing induction of C4 photosynthesis in transitions from dark to high light
The dynamic model presented here extends a previous C4 metabolic model (Wang et al., 2014a,b) to include post-translational regulation, temperature responses of enzymes, dynamic stomatal conductance and leaf energy balance ( Figure 1). The model was built by superimposing the dynamic regulation of enzyme activation and stomatal conductance on the metabolic NADP-ME C4 leaf photosynthesis model of Wang et al. (2014). This was initially parameterized from the literature (Table 2). During induction some C4 metabolic pools, in particular malate in the bundle sheath cytoplasm, rise to very high concentrations (Leegood, 1997). To assess the role of this accumulation of photosynthetic metabolites during induction, the model was first run assuming that all enzymes were fully activated and the stomata open. Termed scenario 1, this resulted in a rapid induction to near steady state within 120 sec (Figure 2a). The major limitation over this period was the time taken for C4 metabolites to accumulate and approach steady state, lagging C3 metabolites ( Figure S1). Leakiness (ϕ), that is the proportion of CO 2 released by decarboxylation in the bundle sheath that diffused back to the mesophyll, reached a minimum at 30 sec, gradually climbing to a steady-state value of approximately 0.22 at approximately 600 sec, indicating that the flux through the C4 cycle continued to limit photosynthesis ( Figure 2b). This limitation was affected by the activity of mutase and enolase, the enzymes that convert PGA to PEP. Increasing the maximum activity of mutase and enolase accelerated induction in scenario 1 ( Figure S2).
In scenario 2, the regulation of PPDK by its regulatory protein (PDRP) substantially slowed the rate of induction (dA/dt) of the CO 2 assimilation rate (A) ( Figure 2a) by limiting the synthesis of phosphoenolpyruvate (PEP), and thus lowering the predicted ϕ ( Figure 2b). In scenario 3, Rubisco regulation alone is added and resulted in a similar decrease in the rate of induction of (dA/dt) to that of scenario 2. It reduced the final steady-state level of A, as a greater proportion of Rubisco now remained inactive (Figure 2a). As would be expected, in contrast to scenario 2, leakiness is high in scenario 3 throughout induction, as the C4 cycle is delivering CO 2 to the bundle sheath, but Rubisco is not fully activated and so is less able to use the CO 2 being released by malate decarboxylation (Figure 2b). Combining the activation of both PPDK and Rubisco to give scenario 4 results in a yet slower rate of induction ( Figure 2a), but the closer coordination of the activation of the two enzymes results in less bundle-sheath leakiness during induction. The simulated leakiness increases for the first 600 sec and then declines to a steady-state value of approximately 0.28 by 1200 sec, reflecting the predicted faster activation of PPDK than Rubisco ( Figure 2b). The addition of dynamic control of the other light-activated enzymes of photosynthetic carbon metabolism ( Figure 1) in scenario 5 produces dynamic responses of A and ϕ that are almost identical to those of scenario 4 ( Figure 2). Finally, superimposing the response of stomatal conductance (g s ) in scenario 6 on the dynamics of A and intercellular CO 2 concentration (C i ) further slows the speed of induction, but dampens the bundle-sheath leakage that would otherwise occur (Figure 2a,b).
The model is shown to predict typical dynamic responses of A and g s , both with respect to pattern and magnitude during induction. The simulation predicts PPDK activation, Rubisco activation and stomatal dynamics as the major limitations, whereas the activation of other enzymes of carbon metabolism and metabolic pool size adjustment had little effect ( Figure 3). The concentration of PDRP regulates the initial phase of the photosynthetic

Measured photosynthetic induction of maize, sorghum and sugarcane
To further analyze the limitations for different C4 crop species, the steady-state and dynamic gas exchange of three major C4 crops were measured, using one widely grown or well-studied genotype for each species: maize B73, sugarcane CP88-1762 and sorghum Tx430. During the transition from dark to high light, the CO 2 assimilation rate of maize rose the fastest, followed by sorghum and then by sugarcane ( Figure 4a). The time required to reach 50% of the steady-state rate (IT50) for the three crops was 196, 237 and 316 sec, respectively. The average CO 2 assimilation rate in the 30-min induction was reduced by 17.7, 20.6 and 24.2% in maize, sorghum and sugarcane, respectively, compared with the steady-state CO 2 assimilation rate. However, maize had a slower rate of increase in g s compared with sorghum and sugarcane ( Figure 4c; Table 2). The C i value dropped rapidly in the first approximately 100 sec, and then slowly increased to the steady-state level. The lowest C i values were 66, 86 and 107 µmol mol −1 in maize, sorghum and sugarcane, respectively, which are 53, 21 and 22% lower than their steady-state C i values (Figure 4b). The low C i values would appear to be insufficient to fully saturate photosynthesis from about 180 to 600 sec after illumination began. Maize showed the highest intrinsic water-use efficiency (iWUE) in the first 600 sec, whereas sorghum had the highest iWUE after 600 sec (Figure 4e). iWUE is the ratio of A to g s . Non-photochemical quenching (NPQ) of the three species rose to a peak at approximately 60 sec and then declined to a steady state at approximately 600 sec, largely in concert with assimilation ( Figure 4d).  (Wang et al., 2014). Here, only enzymes that are light modulated and therefore modified in this new dynamic model are indicated. Green rectangles are driving environmental variables affecting enzyme activities (purple) and stomatal conductance. Blue blocks are state variables calculated from leaf energy balance for T leaf , dynamic stomatal response model for stomatal conductance (g s ), and from the external [CO 2 ], g s and predicted leaf CO 2 uptake rate for C i .  (Wang et al., 2014), which assumes steady-state light activated enzyme activity and stomatal conductance from time zero, that is, throughout. In scenarios 2 and 3, DyPPDK and DyRubisco are added to the steady-state model of scenario 1 to model the induction responses of pyruvate phosphate dikinase (PPDK) and Rubisco, respectively, regulated by the action of the PPDK regulatory protein (PDRP) and by the action of Rubisco activase (Rca) for Rubisco, respectively. Scenario 4 combines scenarios 2 and 3, and scenario 5 includes all light-regulated enzymes. Scenario 6 superimposes stomatal opening on scenario 5. The induction simulates transfer from darkness to full sunlight: 1800 µmol m −2 sec −1 . The input parameters are those listed in

Model parameterization and validation
Using the measured steady-state and dynamic gas exchange data (Figures 4, 5 and Figure S6), we estimated the following photosynthetic parameters: maximum Rubisco activity (V cmax ); maximum PEP carboxylase activity (V pmax ); rate constants for stomatal conductance on darkto-light and light-to-dark transitions (k i and k d , respectively); time constant of Rubisco activation (τ Rubisco ); Figure 3. Simulated induction of leaf CO 2 uptake (A) from dark to high light with variation in (a) the concentration of pyruvate phosphate dikinase regulatory protein (PDRP), (b) the time constant for the activation of Rubisco by Rubisco activase (τ Rubisco ), and (c) the speed of stomatal opening (g s_ki ); k i is the rate constant for increase in stomatal conductance. After dark adaptation, and at time 0 above, the light intensity was raised to 1800 µmol m −2 sec −1 to initiate induction. The input parameters and 'Original values' are listed in Table 2; symbols are defined in Table 1 Table 2). Using these speciesspecific parameters alone, the model was able to closely replicate the measured dynamics of A and g s in all three C4 crops under fluctuating light ( Figure 5). This consisted of 30 min of dark adaptation, followed by 30-min intervals of high light, low light and high light again ( Figure 5).

Factors limiting the speed of photosynthetic induction
Sensitivity analysis of PDRP and Rca concentrations and the speed of the stomatal response indicated that all three limit the rate of photosynthetic induction upon transition from dark to light in the three C4 crops ( Figure 6). However, the strength of each limitation varied between species and with time into induction. In maize B73, sensitivity analysis suggested that PDRP exerts the highest limitation in the first 200 sec of induction, followed by stomatal opening over the next 400 sec and then Rca slightly limits the remaining phase of induction (Figure 6a). In sorghum Tx430 and sugarcane CP88-1762, the concentration of PDRP limits the rate of induction in the first 240 sec, a little longer than in maize (Figure 6b,c).
Rca exerts more influence in sorghum, with a peak at around 420 sec (Figure 6b), whereas the Rca limitation in maize and sugarcane remains approximately constant over this time period (Figure 6a,c). Stomatal limitation was greater in maize and sugarcane compared with sorghum ( Figure 6). In general, PPDK and Rubisco have high control coefficients in the first few minutes; whereas that of PPDK then declines, Rubisco continues to exert control through the mid and final stages of the induction. PEPC also has a high control coefficient from the middle of the induction in sugarcane ( Figure 7e). PPDK and ME have some control in the later stage in maize and sorghum ( Figure 7c). A control coefficient in the present context quantifies the relative influence of a single metabolic step on the rate of CO 2 assimilation. A control coefficient of 1 indicates that the step has complete control and a control coefficient of zero indicates that the step has no control. Except for Rubisco, other light-regulated enzymes of the Calvin-Benson cycle, including glyceraldehyde-3-phosphate dehydrogenase (GAPDH), sedoheptulose-bisphosphatase (SBPase) and phosphoribulokinase (PRK), exerted only mild control in the first 150 sec of the induction (Figure 7b,d,f). . Measured gas exchange parameters upon transition from dark to high light (1800 µmol m −2sec −1 ) after 30 min of dark adaptation: (a) leaf CO 2 uptake rate (A); (b) intercellular CO 2 concentration (C i ); (c) stomata conductance (g s ); (d) nonphotochemical quenching (NPQ); and (e) intrinsic water-use efficiency (iWUE). Bars represent AE1 SEMs for six plants. Leaf CO 2 uptake (A) of the youngest fully expanded leaf was measured on 30day-old maize B73, sugarcane CP88-1762 and sorghum Tx430 plants with a gas exchange system (LI-6800, LI-COR).
© 2021 The Authors. Predicted CO 2 leakiness (ϕ) during photosynthetic induction Predicted ϕ shows an increase as PPDK becomes activated and then declines as simulated Rubisco activity catches up for the three C4 crops. This suggested a loss of coordination between the C4 and Calvin-Benson cycles during induction ( Figure 8). The simulated ϕ of sorghum declines more slowly than that of maize and sugarcane, because of the slower rate of Rubisco activation ( Figure S4; Table 2).

DISCUSSION
The energy-use efficiency of C4 leaves is impacted under fluctuating light In this study, the average photosynthetic rate through the 30-min induction was reduced by 18, 21 and 24% in maize, sorghum and sugarcane, respectively, as compared with the steady-state photosynthetic rate (Figure 4a). This reduction has a very significant effect on the energy-use efficiency and net carbon assimilation of crops in the field, as clouds, wind and the movement of the sun cause frequent light fluctuations within leaf canopies (Kaiser et al., 2018;Tanaka et al., 2019;Wang et al., 2020;Zhu et al., 2004). In sunny conditions, C4 crops have a higher light energy-use efficiency compared with C3 crops as a result of the CO 2 concentrating mechanism, which largely removes the energy cost of photorespiration (Zhu et al., 2010). However, C4 crops may be less resilient to fluctuating light levels, resulting in a decrease in productivity in dynamic light environments (Kubásek et al., 2013). This indicates a significant potential for yield improvement of C4 food and biofuel crops by engineering or breeding for improved speeds of adjustment to fluctuating light. Including post-translational regulation, the temperature response of enzyme activities, dynamic stomatal conductance and a leaf energy balance module, the model closely simulated the measured photosynthetic responses of these crops under fluctuating light ( Figure 5), in contrast to the original metabolic model (Figure 2). This suggests that the model captured the key factors affecting the speed of induction upon light fluctuations ( Figure 5). Using this model, we determined the factors influencing the speed of induction. With the species-specific input parameters (Table 2), the model is able to predict the limiting factors under conditions of fluctuating light (Figures 6 and 7). This has identified potential targets for the improvement of energy-use efficiency in maize, sorghum and sugarcane. Namely, the coordinated upregulation of Rca and the PPDK regulatory protein, as well as increased rates of stomatal adjustment.

Limiting factors during photosynthetic induction
In C3 plants, the rate of photosynthetic induction is mainly limited by the activation of Rubisco, the activation of the enzymes affecting RuBP regeneration and the speed of stomata opening, with the major limitations varying between species (Acevedo-Siaca et al., 2020; De Souza et al., 2020;McAusland et al., 2016;Taylor and Long, 2017). However, the limitations to C4 photosynthetic efficiency under fluctuating light have received little attention.
Through a combination of model simulation and gasexchange experiment, we identified the following limiting factors in photosynthetic induction: (i) the accumulation of C4 photosynthetic intermediates to drive intercellular flux; (ii) the activation of PPDK; (iii) stomata opening; and (iv) the activation of Rubisco.
In our simulations, C4 cycle metabolites took longer to reach a steady state compared with Calvin-Benson cycle enzymes ( Figure S1c,d), although the influence on the induction of photosynthesis was limited to the first 120 sec ( Figure S1a). Also, accelerating the exchange of metabolites between the Calvin-Benson cycle and the C4 cycle, that is, increasing the activity of mutase and enolase, which catalyze the conversion of PGA to PEP, can further reduce the limitation of metabolites during this initial period of induction ( Figure S2). We noted that mutase and enolase exerted higher control at the beginning of the induction and dropped to zero after about 60 sec, based on Figure 6. Simulated changes in the sensitivity coefficients of key parameters through photosynthetic induction for (a) maize, (b) sorghum and (c) sugarcane. After dark adaptation, light intensity was raised to 1800 µmol m −2 sec −1 ; time 0 above. To determine which steps in the system exert the strongest control on dynamic photosynthesis rate, a sensitivity analysis was performed by varying each parameter AE1%. Sensitivity coefficients are calculated as the ratio of change in the value of the parameter divided by change in leaf CO 2 uptake rate (A), individually. If a 1% change in a parameter results in a 1% change in A, the sensitivity coefficient is 1, whereas if the change in A is zero, then the sensitivity coefficient is 0, meaning that no effect is exerted by that parameter. our control analysis for the three C4 crops (Figure 7). However, if the leaves experience a short-term sun fleck, increased rates of photosynthetic metabolite accumulation would improve efficiency. The high concentration of C4 metabolites in NADP-ME species at light saturation results in a slower decline in leaf CO 2 uptake upon transitions from high light to shade, as the decarboxylation of malate continues to provide NADPH, compensating for the decline in NADPH from whole-chain electron transport for a few minutes (Stitt and Zhu, 2014). Our results infer that increasing the concentration of PDRP will increase the photosynthetic efficiency of these C4 plants under the fluctuating light conditions of field crop canopies. This is based on the simulation using the dynamic model developed in this study, which suggested that the concentration of PDRP is a major limitation during the first 180 sec of induction for maize, and for roughly the first 250 sec of induction for sorghum and sugarcane (Figure 6). It regulates both dark-induced inactivation and light-induced activation of PPDK by catalyzing the reversible phosphorylation of a threonine residue (Ashton et al., 1984;Budde et al., 1985;Hatch, 1983, 1985;Chastain, 2010;Chastain et al., 2018). Although these studies elucidated the molecular mechanism for the activation of PPDK by PDRP, the direct effect of PDRP on photosynthetic efficiency has not been analyzed previously. The present analysis suggests that overexpression of PDRP would increase photosynthetic efficiency under field conditions. The sensitivity coefficient of the time constant of stomata opening (k i ) indicated that the speed of stomatal opening was rate limiting from 180 to 600 sec after illumination in maize ( Figure 6). This differed from two previous studies indicating that stomatal conductance was not limiting in maize because C i was always higher than 100 µmol mol −1 during photosynthetic induction (Furbank and Walker, 1985;Usuda and Edwards, 1984). A meta-analysis of responses of A to C i across a number of studies in maize indicated that A was only CO 2 saturated at C i ≥ 100 µmol mol −1 (Pignon and Long, 2020). In the measurements here, C i dropped as low as 66 µmol mol −1 , suggesting that g s was a limitation (Figure 4). However, our study used a higher inducing light intensity of  (Usuda and Edwards, 1984) and 115-1150 µmol m −2sec −1 (Furbank and Walker, 1985), and a longer dark treatment time of 30 min, in comparison with 10 and 20 min, respectively. The higher the light intensity used, the lower C i appeared during the induction (Furbank and Walker, 1985). The longer dark treatment time used here, was to allow sufficient time for stomata to close and for Rubisco to deactivate.
The present analysis indicated that the activation of Rubisco by Rca is the most important limiting factor after the first few minutes of induction, especially in sorghum with slower Rubisco activation (Figures 5 and 6). In Oryza sativa (rice), Rca has been demonstrated to play a crucial role in the regulation of non-steady-state photosynthesis (Yamori et al., 2012). Rubisco is arguably the major limiting enzyme of light-saturated C4 photosynthesis (von Caemmerer, 2000;von Caemmerer et al., 2005;Furbank et al., 1997;Kubien et al., 2003;Wang et al., 2014), and increasing both Rca and Rubisco content have been shown to increase grain yield in maize (Salesse-Smith et al., 2018;Yin et al., 2014). Hence, based on our simulation and previous studies, increasing the activity of Rubisco and Rca in tandem will increase the photosynthetic efficiency under constant and fluctuating light.
Phosphoenolpyruvate carboxylase (PEPC) appears not to restrict photosynthesis under steady-state conditions, except under conditions inducing a low C i , such as drought (Pignon and Long, 2020). However, as C i dropped below 100 µmol mol −1 during induction (Figure 4b), sensitivity analysis indicated that increasing PEPC would increase photosynthetic efficiency from 180 sec to about 600 sec during induction in maize and sorghum (Figure 7a,c). In sugarcane, however, PEPC limits the steady-state photosynthetic rate of sugarcane, through its lower V pmax compared with maize and sorghum ( Figure 7e; Table 2). Furbank et al. (1997) concluded from antisense manipulations that PPDK and Rubisco shared metabolic control of steady-state light-saturated photosynthesis in the C4 dicot Flaveria bidentis. The limited studies of C4 photosynthesis under fluctuating light have focused on maize. Two early studies indicated that photosynthesis reached a maximum rate after about 15-25 min in maize (Furbank and Walker, 1985;Usuda and Edwards, 1984), which is comparable with the observations and simulations here (Figures 4a and 5).

Differences in the limiting factors of photosynthetic induction among species
This study is limited to single accessions of each of three NADP-ME C4 species. Therefore, it cannot be generalized to the species studied here. However, the examination of individuals from three distinct species of the monophyletic Andropogonae, all C4-NADPME plants, is likely to have revealed limitations that apply across this key clade of food and energy crops. They therefore point to manipulations that could improve photosynthetic efficiency and yields across the clade. Although there were many similarities, some differences were found. Maize, as perhaps the species most intensively bred for productivity, showed the fastest induction and greatest efficiency of carbon gain over the period of induction, whereas sugarcane was the slowest (Figure 4). Whether these are species characteristics can only be determined by analyzing a wider range of genotypes of each crop. Characterizing within-species variation would also show the potential for improving non-steady-state photosynthesis through breeding. In rice, intraspecific genetic variation in non-steady-state photosynthetic efficiency was found to be substantially greater than in the steady state, suggesting an overlooked target for improvement that might similarly be available in these crops (Acevedo-Siaca et al., 2020a,b).
Maize showed the fastest induction, because of more PDRP and faster τ Rubisco (Table 2), which indicated that maize has faster PPDK and Rubisco activation capacity. However, the stomatal response of maize is slow (Table 2). Here the stomata were one of the major limiting factors during the induction process ( Figure 6). This conclusion is different from previous studies (Furbank and Walker, 1985;Usuda and Edwards, 1984), and the possible reasons were explained in the previous section. Speeding up stomatal opening and closing is the key to speed up the photosynthetic response while maintaining water-use efficiency. New combined thermal and modulated fluorescence techniques now provide a potential high-throughput means to screen germplasm for this trait (Pignon et al., 2021;Vialet-Chabrand and Lawson, 2019;Vialet-Chabrand and Lawson, 2020). Bioengineering for more and smaller stomatal complexes would be another approach (Drake et al., 2013). For sorghum, the speed of stomatal opening had little effect on A during induction ( Figure 6). Enzyme activities were the main limiting factors, that is, the concentration of PDRP (i.e. the activation rate of PPDK), the activation rate of Rubisco (τ Rubisco ) and Rubisco activity (V cmax ) (Figures 6  and 7). Thus, increasing the activity of PDRP, Rca and Rubisco would lead to higher dynamic photosynthesis. However, analysis of water-use efficiency across a wide range of sorghum germplasm suggests at the species level that the speed of stomatal adjustment is also important (Pignon et al., 2021).
For sugarcane, its dynamic photosynthetic efficiency was co-limited by many factors, including the rate of stomatal opening, the concentration of PDRP, the activation rate of Rubisco (τ Rubisco ) and Rubisco activity (V cmax ). In addition, a high control coefficient of PEPC, relative to the other species, was found in sugarcane during induction and in the steady state ( Figure 7). Therefore, to improve dynamic photosynthesis, all the limiting factors above may need to be considered. However, with only one accession of each species, our study cannot determine whether these are actual species differences or simply the result of the accessions that were chosen.

Imbalances in the regulation of C3 and C4 cycles
Coordination between the C3 and C4 cycles is essential to the efficiency of C4 photosynthesis. Leakiness (ϕ) describes the proportion of carbon fixed by PEPC that retrodiffuses back out of bundle-sheath cells into the mesophyll (Equation 30). It was estimated to be about 0.2 in several C4 species under various environmental conditions (Henderson et al., 1992) and between 0.20 and 0.22 in a recent study of maize (Salesse-Smith et al. 2018). Our simulated steadystate ϕ values of the three species lie between 0.2 and 0.3 ( Figure 8). In our simulation, the leakiness is predicted to change during the induction through an imbalance in the regulation of the Calvin-Benson cycle and the C4 cycle, especially when the activation of Rubisco is slower (Figure 8, Sorghum). Thus, increasing the activation speed of Rubisco is the first choice to improve efficiency during transitions from dark to light. Although the activation of PPDK is also one of the main limiting factors, increasing the PDRP concentration could not significantly improve photosynthetic efficiency without a larger increase in the speed of Rubisco activation. Overall, this study has identified several potential opportunities for increasing photosynthetic efficiency in these major crops during the frequent light fluctuations that occur in field canopies.

Model development
We developed a generic dynamic systems model of C4 photosynthesis from the previous NADP-ME metabolic model for maize (Wang et al., 2014a,b, Appendix S2). The NADP-ME metabolic model is an ordinary differential equation model including all individual steps in C4 photosynthetic carbon metabolism. Here, we extended the model to include the post-translational regulation and temperature response of enzyme activities, together with the dynamics of stomatal conductance and leaf energy balance. The model was implemented in MATLAB. To save space and make the text more readable, abbreviations have been introduced throughout this section. They can be seen in Table 1.
Post-translational regulation of enzyme activity. PPDK activation state-The activity of PPDK is regulated by PDRP, which is affected by the level of incident light via ADP concentration (Ashton et al., 1984;Burnell and Hatch, 1985;Chastain, 2010). PDRP is a bifunctional protein kinase/protein phosphatase, catalyzing the reversible phosphorylation of PPDK. The inactivation rate (V PDRP I Þ and activation rate (V PDRP A Þ were calculated by the following equations: where [PDRP] Mchl is the PDRP concentration in the mesophyll cell chloroplasts, and k cat PDRP I and k cat PDRP A are the turnover number of PDRP for the inactivation and activation reaction, respectively. E ½ Mchl is the concentration of active PPDK in the mesophyll chloroplasts; EP ½ Mchl is the concentration of inactive PPDK in the mesophyll chloroplasts.
Differential equations, parameters and their sources are listed in Appendix S1.
Rubisco activation state-The time constant of Rubisco activation was determined from the measured kinetics of photosynthetic gas exchange following transitions from dark to high light, using the method given by Woodrow and Mott (1989) where τ Rubisco is the rate constant of Rubisco activation catalyzed by Rca. V max_Rubisco_i is the transient maximal Rubisco activity; V max_Rubisco_s is the steady-state maximal Rubisco activity, which is related to the Rca concentration, [Rca] (Mott and Woodrow, 2000). The total [Rca] is calculated using measured τ Rubisco (Equation 25; Table 2): where k is a constant, which is 216.9 min mg m −2 (Mott and Woodrow, 2000). Steady-state maximal Rubisco activity is calculated with the following equations: where V max Rubisco is the theoretical maximum activity of Rubisco.
[Rca] A is the concentration of active Rca, which is regulated by light intensity (Appendix S1). K activase is a constant, which equals 12.3 mg m −2 (Mott and Woodrow, 2000).

Activation of enzymes regulated via light intensity-The
model used a simplified equation for the light regulation of ATP synthase (ATPase), fructose-1:6-bisphosphatase (FBPase), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), phosphoribulose kinase (PRK), rubisco activase (Rca) and sedoheptulose-1:7-bisphosphatase (SBPase): where V max_E_i is the transient maximal enzyme activity, τ E is the rate constant of the activation of each enzyme and V max E s is the steady-state maximal enzyme activity, as affected by light intensity (I). k E_A and c E_A are two constants, that is, the slope and intercept of the linear relationship of the proportion of activated enzyme (a E s ) as a function of I. V max E is the activity of the enzyme when fully activated. Although the activation of PEPC is regulated by light via phosphorylation, the whole pathway and parameters of this regulation have not been measured quantitatively. Thus, the dynamics of  Temperature response of enzymes. In order to simulate the effects of fluctuating leaf temperature with fluctuations in light, the Arrhenius equation (Johnson et al., 1942) and Q 10 function were used to adjust the enzymatic parameters to the actual leaf temperature (T leaf ). The formula used for each parameter was determined based on the availability of experimental data. The temperature response of the maximum activity of carbonic anhydrase (CA) and PEPC (V max_CA and V max_PEPC ) were incorporated into the model using a peaked Arrhenius function (Johnson, Eyring and Williams, 1942): where E a is the exponential rate of the rise, H d describes the rate of decrease at supraoptimal temperatures and ΔS is the entropy factor. The remperature response of enzymatic parameters of pyruvate phosphate dikinase (V max_PPDK ), electron transport capacity (J max ) and Rubisco (V max_Rubisco_CO2 , V max_Rubisco_O2 /V max_Rubisco_CO2 , K o and K c ) were incorporated into the model using an Arrhenius function: Parameters and sources are listed in Appendix S1.
For other enzymes, a Q 10 function was used to estimate the temperature response of the maximum activity, as described previously (Woodrow and Berry, 1988). Q 10 was set as 2: 10 Enz : Dynamic stomatal response. Ball-Berry model parameters for predicting steady-state stomatal conductance (Ball et al., 1987) were obtained from light response curves measured for each C4 crop evaluated in this study. In the Ball-Berry model, stomatal conductance was with a function of A, relative humidity (RH) and CO 2 concentration at the leaf surface (C a ): where Slope BB is the slope of the relationship between g s_steady and A × RH/C a and Intercept BB is the residual stomatal conductance. Slope BB and Intercept BB were estimated by linear regression of AÁRH Ca and g s_steady from the light response curve (A-Q curve) measurement.
Dynamic stomatal conductance (g s ) was estimated with the following equation: where g s_steady is the steady-state stomatal conductance calculated by the Ball-Berry model (Equation 13;Ball et al., 1987); k (k i or k d ) is the rate constant of the stomata conductance response calculated from the measured stomata dynamics of the three C4 crops, and k i and k d represent the rate constants of stomata conductance increasing and decreasing, respectively ( Dynamic leaf energy balance. For leaf energy balance, the equations used in our model were based on the method of Nikolov et al. (1995). According to this model, leaf energy balance takes account of intercepted short-and long-wave radiation, radiative energy loss from the leaf, convection and latent heat loss in transpiration. The net photosynthesis rate (A), stomatal conductance and leaf temperature are interdependent. For example, A affects stomatal conductance, stomatal conductance affects leaf temperature and leaf temperature affects A. Instead of solving these steady-state circular connections iteratively (Nikolov et al., 1995), a differential equation describes leaf temperature (T leaf ) change (Equation 15): where PAR abs is the absorbed photosynthetic active radiation, assuming that 85% of PAR is absorbed by the leaf, NIR is the absorbed near-infrared radiation and LR is the absorbed longwave radiation. Both NIR and LR were set to zero. C p is the specific heat capacity of the leaf, and here we assumed it is the same as the specific heat capacity of water (4.184 J g −1°C−1 ). m leaf is the specific leaf fresh weight (g m −2 ) and was set as 198 g m −2 for all species based on the measured value of maize leaves (197.9 AE 4.5 g m −2 ). Humidity in the leaf internal air space is assumed to be saturated at the temperature of the leaf. H and LE are the sensible and latent heat fluxes, respectively. E is the emitted long wave radiation and M e is the energy consumed in photosynthesis (Nikolov et al., 1995). The boundary layer conductance to heat is calculated as g bh = 0.924g b (Nikolov et al., 1995). C P_air is the specific heat capacity of air (29.3 J mol −1°C−1 ), C lv is the latent heat of vaporization of water (44 000 J mol −1 ) and g l is the total conductance of the stomata and the boundary layer. ε is the leaf emisivity of long-wave radiation and σ is the Boltzman constant.
Boundary layer conductance. Boundary layer conductance was calculated following Nikolov et al. (1995), and both free and forced convection was considered in determining the boundary layer conductance of the leaf. The leaf boundary layer conductance to vapor transport is the maximum of g bf and g br : The forced-convective and free-convective leaf boundary layer conductance is computed as: g br ¼ c e T 0:56 leafk T leafK þ 120 P a 0:5 ΔT ð Þ 0:25 , where d o is the characteristic dimension of a leaf (leaf width), ΔT is the temperature difference between the leaf and the local air (Monteith and Unsworth, 1990), u is the wind velocity, and c f and c e are two constants. Plant material and growth conditions. Maize B73, sugarcane CP88-1762 and sorghum Tx430 were grown in a controlledenvironment glasshouse at 28°C (day)/24°C (night). Maize and sorghum were grown from seed and sugarcane CP88-176 was grown from stem cuttings. The positions of the plants in the glasshouse were re-randomized every week to avoid the influence of environmental variations within the glasshouse. From 25 July to 8 August 2019, six biological replicates were measured in a randomized experimental design for each species in each measurement.
Steady-state gas exchange measurements and parameter estimation. Leaf gas exchange of the youngest fully expanded leaf was measured on plants at 30-35 days old with a gas exchange system (LI-6800; LI-COR, Lincoln, NE, USA). The leaf chamber temperature was set at 28°C, with a water vapor pressure deficit of 1.32 KPa and a flow rate of 500 µmol sec −1 for all gas exchange measurements. For the response of A to intracellular CO 2 concentration curves (A-C i curves), the leaf was acclimated to a light intensity of 1800 µmol m −2 sec −1 and a CO 2 concentration of 400 µmol mol −1 . After both A and g s reached steady state, the CO 2 concentration of the influent gas was varied in the following sequence: 400,300,200,120,70,40,20,10,400,400,400, 600, 800, 1200 and 1500 µmol mol −1 .
The maximum Rubisco activity (V cmax ) and maximum PEPC activity (V pmax ) were estimated from the A-C i curves using the equations of von Caemmerer (2000). In order to obtain the relationship between estimated V pmax and the theoretical maximal PEPC activity (V max_PEPC ) in our model, and similarly the relationship between V cmax and the theoretical maximal Rubisco activity (V max_Rubisco ), we introduced two variables (f vpmax and f vcmax ) into the simulation: Both f vpmax and f vcmax were estimated by minimizing the sum (S fvPEPC andS fvcmax ) of the squared differences between the estimated A (A e_Ci ) from the dynamic model and the measured A (A m_Ci ), in response to intercellular CO 2 (A-C i curve), using the least-squares method for each species: f vPEPC was estimated using the initial slope (s Am Ci ) of the measured A-C i curve (CO 2 air = 120, 70, 40, 20, 10 µmol mol −1 ); f vcmax was estimated using CO 2 -saturated A m_Ci (CO 2 air = 800, 1200 and 1500 µmol mol −1 ) ( Figure S3). The steady-state V max of the other enzymes of C4 and C3 metabolism in Figure 1 were scaled for each species with f vpmax and f vcmax , respectively.
To define the response of A to light intensity (A-Q curves), the leaf was acclimated to a light intensity of 1800 µmol m −2 sec −1 and a CO 2 concentration of 400 µmol mol −1 . After leaf gas exchange reached the steady state, the light intensity in the chamber was changed in the following sequence: 2000, 1500, 1000, 500, 300, 200, 100 and 50 µmol m −2 sec −1 . The gas exchange data were logged after 5 min to ensure that there was enough time for the transpiration, and therefore stomatal condcutance, to reach the steady state at each light level. Ball-Berry model parameters (Ball et al., 1987) were estimated by the linear regression of AÁRH C a and g s_steady from data from the A-Q curves, including the prediction of steady-state stomatal conductance (g s_steady ) for each species (Equation 13).
Dynamic gas exchange measurements and parameter estimation. Gas exchange during photosynthetic induction was measured in the transition from darkness to high light (1800 μmol m −2 sec −1 ) to determine the kinetics of Rubisco activation in these C4 crops (τ Rubisco ). The leaf was first acclimated to darkness for 30 min, with a CO 2 concentration of 400 µmol mol −1 and then the light intensity was changed to 1800 μmol m −2 sec −1 for 30 min, which was more than sufficient time for leaf CO 2 uptake and stomatal conductance to reach the steady state. Leaf gas exchange was logged before the light was turned on, and then logged every 10 sec for the following 30 min. The time constant of Rubisco activation (τ Rubisco ) was estimated from the linear portion of the semi-logarithmic plot of photosynthesis with time Mott, 1989, 1993; Figure S4). The slope of this portion was determined by the linear regression of the data between 3 and 7 min. The value of τ Rubisco was calculated as: The calculated values of the three C4 species are listed in Table 2.
To further evaluate the response of gas exchange in C4 plants under fluctuating light, following this 30 min of induction the responses to the transition from high to low and back to high light (i.e. relaxation curves followed by induction curves) were measured. This involved decreasing light to 200 µmol m −2 sec −1 PPFD for 30 min and then returning to 1800 µmol m −2 sec −1 PPFD for an additional 30 min. Gas exchange was recorded every 10 sec.
Rate constants were calculated for g s increase on transfer from low light (200 µmol m −2 sec −1 PPFD) to high light (k i ), and again for the decrease in g s on return to low light (k d ). The measured time series for stomatal conductance changes were fitted with the following equation (Vialet-Chabrand et al., 2017): where g max is the maximum stomata conductance, g 0 is the minimum stomata conductance, t is time and k (k i or k d ) is the rate constant of g s . g max , g 0 and k were estimated using Equation 28 by the curve fitting function (FIT) in MATLAB (MathWorks, https:// www.mathworks.com). Mitochondrial respiration (R d ) was estimated from the measured CO 2 efflux after 30 min of dark adaptation. The concentration of PDRP was estimated by minimizing the difference between the estimated A (A e_t ) from the dynamic model and the measured A (A m_t ) at the beginning of the induction using the least-squares method, which minimizes the sum (S PDRP ) of the squared difference between estimated and measured A in the beginning of the photosynthetic induction ( Figure S5). Data points from 1-3 min of the induction were used.
Model parameterization. The model took the following 11 photosynthetic parameters estimated from measured gas exchange data as input variables: maximum Rubisco activity (V cmax and f vcmax ); maximum PEPC activity (V pmax and f vpmax ); the rate constant of stomata conductance increase and decrease (k i and k d ); the time constant of Rubisco activation (τ Rubisco ); mitochondrial respiration (R d ); [PDRP]; and the Ball-Berry slope (Slope BB ) and intercept (Intercept BB ) ( Table 2). The estimation methods of the input variables were described in the previous section (Gas exchange measurement and parameter estimation). Parameters and equations affected by the input variables are listed in Appendix S1.

Model prediction
CO 2 uptake rate (A) and leakiness (ϕ) calculation. During the simulation, metabolite concentrations and reaction rates were extracted from the model. The velocity of CO 2 flowing into the leaf via stomata was used to represent A. Leakiness (ϕ) describes the proportion of carbon fixed by PEPC that subsequently leaks out of the bundle sheath cells. Thus, ϕ was calculated as: where the CO 2 leak rate (v CO2_leak ) is determined by the permeability of CO 2 through plasmodesmata (P CO2_pd ) and the concentration gradient of CO 2 between bundle sheath cytosol and mesophyll cytosol, [CO 2 ] BSC -[CO 2 ] MC , and v PEPC is the velocity of carbon fixation by PEPC.
Sensitivity analysis. The sensitivity coefficient (SC p ) gives the relative fractional change in the simulated result with a fractional change in the input variable (p), SC p is the partial derivative used to describe how the output estimate varies with changes in the values of the input parameter (p), where the output in this study is the estimated leaf CO 2 uptake rate (A): where the variable (p) was both increased and decreased by 1% individually in the model to calculate the new value of A (A + and A − , respectively) in order to identify the parameters influencing A.
The flux control coefficient of each enzyme (FCC) was also estimated by Equation 31, using the maximal activity (V max_E ) of the enzyme as the variable (p). Figure S1. Simulated photosynthetic induction using metabolic model without post-translational regulation of enzymes and delay of stomata conductance. Figure S2. Estimated influence of mutase and enolase on photosynthetic induction. Figure S3. Estimation of f vPEPC and f vRubisco using least-squares method. Figure S4. Semilogarithmic plot of the difference between photosynthesis (A) and maximum photosynthesis (A f ) as a function of time. Figure S5. Estimation of PPDK regulatory protein concentration, [PDRP], using measured photosynthetic induction curves. Figure S6. Measured CO 2 response curves and light response curves of maize B73, sorghum Tx430 and sugarcane CP88-1762. Figure S7. Calculated Ball-Berry slope and intercept using gas exchange data from the light response curves of maize, sorghum and sugarcane.