Can increased leaf photosynthesis be converted into higher crop mass production? A simulation study for rice using the crop model GECROS

Highlight We model the impact of improving photosynthesis on rice productivity. We consider all major photosynthesis-enhancing approaches under one umbrella and scale them from leaf photosynthesis to crop production.


Introduction
Yields of major crops have increased steadily during the last decades. Fischer et al., (2014) claimed that, in order to ensure food and energy security for a growing and increasingly demanding population, staple crop production will need to grow by 60% from 2010 to 2050, with the greatest increases in the next 20 years. They further stressed that higher rates of increase than the current rate should be aimed for to drive faster reductions in world hunger and to guard against unanticipated negative contingencies, for example as a result of increasing frequencies of extreme weather under global climate change.
Crop yield per area of land is the production of mass per unit area multiplied by harvest index. Yield gains associated with the first Green Revolution in cereal crops such as wheat and rice were mainly due to increased harvest index by introducing (semi-)dwarfing genes (e.g. Miflin, 2000;Sadras and Lawson, 2011). For further progress to be made, improvement in crop mass production via increasing leaf and canopy photosynthetic capacity and efficiency should be explored Murchie et al., 2009;Parry et al., 2011;Ort et al., 2015). Evidence suggests that genetic variation in leaf photosynthesis has not been exploited to be incorporated into crop cultivars (e.g. Driever et al., 2014), except for recent releases in which yield gains were accompanied, to some extent, by traits related to increased leaf photosynthesis (Fischer et al., 2014).
Arguably, exploiting natural genetic variation is still the most feasible approach to improve yield traits including photosynthesis (Flood et al., 2011). For example, exploring natural variation in mesophyll conductance for CO 2 diffusion (g m ) may improve photosynthesis (Gu et al., 2012;Flexas et al., 2013;Chen et al., 2014). However, very often there is little correlation between leaf photosynthesis and crop productivity across germplasm (Driever et al., 2014;Koester et al., 2016) or across individual lines of a segregating population (e.g. Gu et al., 2014a, b), partly because natural variation in leaf photosynthesis and underlying traits is generally small (Driever et al., 2014). To enhance crop productivity at a greater pace, genetic engineering and synthetic biology approaches to improving leaf photosynthesis should be explored (Long et al., , 2015Singh et al., 2014;Ort et al., 2015;. Major crops such as rice follow the pathway of C 3 photosynthesis. Compared with C 4 crops such as maize, C 3 crops have lower photosynthetic productivity primarily because ~20-35% of the carbohydrate is lost through photorespiration Walker et al., 2016), resulting from the oxygenation of ribulose-1,5-biphosphate by Rubisco, the primary enzyme for CO 2 fixation. Various genetic engineering routes to enhance C 3 leaf photosynthesis have been proposed to suppress photorespiration, thereby improving crop productivity. These include: replacing Rubisco with foreign forms with higher specificity for CO 2 relative to O 2 (S c/o ) (Whitney et al., 2001;Zhu et al., 2004;Parry et al., 2011), designing a photorespiratory bypass (Kebeish et al., 2007), and transforming the C 4 CO 2 -concentrating mechanism (CCM) (e.g. von Caemmerer et al., 2012) or the cyanobacterial bicarbonate-based CCM (Price et al., , 2013Lin et al., 2014) into main C 3 crops.
Some of the engineering approaches have already made progress and were evaluated experimentally for mass production using model plants, such as Arabidopsis thaliana engineered with a photorespiratory bypass (Kebeish et al., 2007) and tobacco with a cyanobacterial bicarbonate transporter (Pengelly et al., 2014) and with accelerated recovery from photoprotection . Any progress for major crops may not be expected in the near future, and modelling should be considered as an important tool to assess the potential of yield improvement by these photosynthesis-enhancing routes. Many researchers (e.g. Zhu et al., 2004;Song et al., 2013; have published modelling studies to assess the potential benefit of using various routes in improving photosynthesis; but most of their analyses were based on the simulation of canopy photosynthesis at a fixed leaf area index (LAI) for a given environmental condition.
During the growth cycle of annual field crops, LAI expands initially, reaches its maximum size, and then senesces, and complex interactions and feedback mechanisms can occur between photosynthesis and other physiological components (Yin and Struik, 2008). These complexities should be considered when scaling up from instantaneous leaf assimilation to daily canopy photosynthesis and to total mass production over the growing season (Boote et al., 2013). To that end, Gu et al., (2014a) ran numerical simulations using a full crop growth model, GECROS (Yin and van Laar, 2005), and examined the potential of exploiting the natural genetic variation in leaf photosynthesis within a single segregating population for contributing to crop productivity in rice (Oryza sativa L.) under field conditions.
Here we ran the crop model GECROS to quantify the extent to which improved leaf photosynthesis, predominantly from genetic engineering to suppress photorespiration, can result in an expected increase in crop mass production under well-watered, as well as water-limited, field conditions, using rice as an example.

Materials and methods
Model algorithms and approach in applying GECROS (v4.0) for this study are outlined below. Model parameters, if not defined in the text, are given in Table 1.

C 3 photosynthesis model
The model of Farquhar et al., (1980; the FvCB model hereafter) calculates net CO 2 assimilation rate (A) as the minimum of the Rubisco-limited (A c ) and e − transport-limited (A j ) rates. The two limiting rates can be expressed collectively as: where for A c , x 1 = V cmax and x 2 = K mC (1 + O/K mO ); for A j , x 1 = [1f pseudo /(1-f cyc )]J 2 /4 and x 2 = 2Oγ * , where x 1 is written according to the FvCB model extended by Yin et al. (2004) to be compatible with a C 4 model for which accounting for f cyc is required (see later). In the model, C c and O are the CO 2 and O 2 level, respectively, at the carboxylation sites of Rubisco, J 2 is the total PSII e − transport rate, and Γ * , defined as Oγ * (where γ * is half of the inverse of S c/o ) is the CO 2 compensation point in the absence of day respiration (R d ).
The submodel for stomatal conductance for CO 2 transfer (g s ) is: where g 0 is the residual value of g s when irradiance approaches zero, C i* is the intercellular CO 2 level (C i ) at which A + R d = 0, and f vpd is the relative effect of leaf-to-air vapour difference (VPD) on g s (see later). CO 2 transfer from C a (the ambient CO 2 level) to C c can be written as (Flexas et al., 2013):  (Cousins et al., 2010;Perdomo et al., 2015) χ Jmax25 Linear slope of maximum PSII e − transport rate at 25 °C (J max25 ) versus (n-n b ) (μmol s −1 g −1 ) 100 Harley et al., (1992);   0.01 Leuning (1995) 0.01 Assumed to be the same as for C 3 a 1 Empirical constant for g s response to VPD (-) 0.9 Derived from Morison and Gifford (1983) 0.9 Set the same as for C 3 crops c b 1 Empirical constant for g s response to VPD (kPa −1 ) 0.15 Derived from Morison and Gifford (1983) 0.15 Set the same as for C 3 crops c χ gm25 Linear slope of mesophyll conductance at 25 °C (g m25 ) versus (n-n b ) (mol s −1 g −1 ) 0.125 Derived from data of   Table 2). b Where n is leaf nitrogen (g N m −2 ); and n b is the base leaf nitrogen, below which no leaf photosynthesis is observed. c Data of Morison and Gifford (1983) showed that stomatal sensitivity to VPD could differ between C 3 and C 4 ; such a difference can be mimicked by our stomatal conductance model, Equation 2 for C 3 and Equation 11 for C 4 leaves, when using the same values of a 1 and b 1 . d Parameter set in GECROS to be dependent on crop species; the value 88 380 was set as default for rice (Yin and van Laar, 2005). Leaf photosynthesis and rice productivity | 2349 Combining Equations 1-4 gives a standard cubic equation, as shown in Supplementary Text 1 at JXB online.

C 4 photosynthesis model
The C 4 model of von Caemmerer and Furbank (1999), as modified by Struik (2009, 2012), is used here. In C 4 plants, CO 2 is fixed initially in the mesophyll by phosphoenolpyruvate (PEP) carboxylase into C 4 acids that are then decarboxylated to supply CO 2 to Rubisco, which is localized in the bundle-sheath chloroplasts. The co-ordinated functioning of the 'Kranz' anatomy and C 4 biochemistry enables an effective CCM. The extra ATP consumption for sustaining the CCM requires a higher f cyc in C 4 than in C 3 photosynthesis Nakamura et al., 2013). The rate of PEP carboxylation (V p ) could be limited either by the PEP carboxylase or by the rate of e − transport ): where ε p is the initial carboxylation efficiency of the PEP carboxylase, φ is the extra ATP required for the CCM per CO 2 fixed, and z is the conversion factor of J 2 into the ATP production rate: (here h is the H + :ATP ratio; Yin et al., 2004;Yin and Struik, 2012), and x represents the fraction of ATP partitioned to the reactions associated with the operation of V p . In the standard C 4 model for malic-enzyme subtypes such as crop plants maize and sorghum, x was set to 0.4, arising from φ/(3 + φ), where φ = 2, and 3 is mol ATP required for the Calvin cycle to fix 1 mol CO 2 . An effective CCM requires a small bundle-sheath conductance (g bs ) as g bs determines the CO 2 leakage from the bundle sheath to the mesophyll (L) that affects CO 2 assimilation (von Caemmerer and Furbank, 1999): Equations 5-7 can be combined to result in: ). The rate of CO 2 fixation by Rubisco is modelled in the same way as for C 3 photosynthesis: where x 1 = V cmax , x 2 = K mC /K mO , x 3 = K mC for the enzyme (Rubisco)limited rate, and x 1 = [1-f pseudo /(1-f cyc )]J 2 /4, x 2 = 2γ * , and x 3 = 0 for the e − transport-limited rate. This form of the e − transport-limited rate implies that it is the NADPH supply that causes the e − transport limitation in C 4 photosynthesis, in comparison with the standard C 4 model in which the ATP supply was assumed to cause the e − transport limitation (von Caemmerer and Furbank, 1999). Yin and Struik (2012) discussed the rationale that either the ATP-or the NADPH-limited form can be used for modelling C 4 photosynthesis provided that f cyc and f pseudo are set as appropriate. We prefer to use the NADPH-limited form here because the ATP-limited form gives (Yin et al., 2011), which would predict a monotonic increase in ATP production rate, thus in an e − transport-limited carboxylation rate, with increasing f cyc . This does not agree with the more efficient CCM in terms of ATP use (e.g. cyanobacterial CCM; Price et al., 2011). Using the NADPH-limited form allows a revised C 4 model to simulate photosynthesis of other CCM systems (see below) and to be consistent with the C 3 photosynthesis modelling where the NADPH-limited form is predominantly used. A relationship for O 2 partial pressure between the intercellular air space (O i ) and the sites around Rubisco in bundle-sheath cells (O) is described as (von Caemmerer and Furbank, 1999): A model for g s of C 4 leaves was formulated in a way that slightly differed from Equation 2 of the C 3 counterpart, to solve analytically for A in C 4 photosynthesis (Yin and Struik, 2009): All C 3 default parameters in Table 1 but χ gm25 = 0.375 2 Improved Rubisco specificity S c/o a C 3 All C 3 default parameters in Table 1 but S c/o = 4427 3 Improved value for both g m and S c/o C 3 All C 3 default parameters in Table 1 but χ gm25 = 0.375 and S c/o = 4427 4 C 4 biochemistry introduced C 4 All C 4 parameters (including χ Vcmax25 and χ Jmax25 b ) in Table 1, but χ bs25 = 0.125 5 C 4 Kranz anatomy introduced effectively to minimize CO 2 leakage C 4 All default C 3 enzymatic parameters plus necessary C 4 parameters to run C 4 model in Table 1, but χ bs25 = 0.007 6 Complete C 4 mechanism engineered C 4 All C 4 parameters in Table 1, including low χ bs25 (= 0.007) 7 Only cyanobacterial bicarbonate transporters engineered C 4 All C 3 default parameters plus necessary C 4 parameters to run C 4 model in Table 1, but χ gbs25 = 0.125, φ = 0.75, x = 0.2, f cyc = 0.18, and R m = R d 8 More elaborate cyanobacterial CCM added C 4 The same as Route 7, but χ gbs25 = 0.007 9 Complete cyanobacterial CCM engineered C 4 The same as Route 8, but with χ Vcmax25 = 93 and χ Jmax25 = 200 c a This route assumes that crop plants are engineered to have a high S c/o25 of the non-green alga Griffithsia monilis while maintaining a similar Rubisco turnover rate (Whitney et al., 2001); any effect of the trade-off between Rubisco S c/o and carboxylase turnover rate was not quantified here, and readers are suggested to refer to Zhu et al., (2014) on this effect. b Based on measurements on existing maize and wheat plants, parameters χ Vcmax25 and χ Jmax25 have higher values in C 4 than in C 3 leaves (Table 1), probably reflecting the acclimation of C 4 enzymatic activities to high a CO 2 environment within the bundle-sheath compartment. While strictly speaking these higher values cannot be guaranteed for hypothetical C 4 plants of Route 4 which is not yet incorporated with the full Kranz anatomy, high values of χ Vcmax25 and χ Jmax25 for maize plants (Table 1) were used here for simulation of Route 4 in order to represent the full package of the C 4 biochemistry components. c Cyanobacterial Rubisco has a higher carboxylation rate than C 3 Rubisco (Hanson et al., 2016), allowing a higher investment of nitrogen in other photosynthetic protein components. However, we are not aware of the N cost for e − transport protein components in cyanobacteria for estimating χ Jmax25 . For simplicity, χ Vcmax25 and χ Jmax25 for maize plants (Table 1) are used for this route, based on the expectation of engineering cyanobacterial CCM that approaches typical C 4 photosynthetic capacities (Price et al., 2013).
where C s is the CO 2 level at leaf surface, and C s* is the C s -based CO 2 compensation point in the absence of R d and can be calculated as [g bs γ * O i -(1 + γ * α/u oc )R d + R m ]/(g bs + ε p ) ).
Equation 11 for C 4 and Equation 2 for C 3 , although both empirical, can reproduce experimentally observed linear relationships between A and g s across various levels of irradiance and nutrients (e.g. Wong et al., 1985) (see Supplementary Fig. S1). Equation 3 also applies to C 4 photosynthesis. Combining Equations 3 and 8-11 can yield the standard cubic equation that gives the prediction of A (Supplementary Text 1).
Algorithms common to C 3 and C 4 photosynthesis Some common algorithms were used for C 3 and C 4 models. First, J 2 is described as a function of absorbed irradiance I abs as (Yin et al., 2006;Yin and Struik, 2012): differs from the equation used in the standard FvCB model, in that f cyc , Φ 2LL , and r 2/1 are introduced. We consider Equation 12 as a better choice as it accounts for the decrease of the overall noncyclic e − transport efficiency (α 2LL ) with increasing cyclic e − transport, which runs at a higher rate in cases involving the CCM.
Secondly, in the g s model, f vpd is the function for the effect of VPD, which may be described phenomenologically as where a 1 and b 1 represent the C i :C a ratio in water vapour-saturated air and the slope of the decrease of this ratio with increasing VPD, respectively, if g 0 in Equation 2 or 11 approaches nil. Thirdly, a number of parameters are related to leaf temperature (T l ), and some of these can be described by the Arrhenius equation normalized with respect to 25 °C: where R is the universal gas constant (8.314 J K −1 mol −1 ). Equation 14 applies to R d , γ * , V cmax , K mC , K mO , and u oc . The temperature response of J max , ε p , g m , and g bs is described by the modified Arrhenius equation: Fourthly, the values at 25 °C of parameters V cmax , J max , ε p , g m , and g bs can be further quantified as a linear function of leaf nitrogen (N) content (n) above a certain base value (n b ): where χ has different values for different parameters (e.g. Harley et al., 1992;Yin et al., 2011). We estimated χ Vcmax25 for C 3 leaves from existing data and then projected to C 4 leaves (Table 1), based on the reported higher catalytic turnover rate of C 4 Rubisco than that of C 3 Rubisco (Seemann et al., 1984;Sage, 2002;Cousins et al., 2010;Perdomo et al., 2015). There is less information about the difference in χ Jmax25 between C 3 and C 4 types, but our χ Jmax25 estimates (Table 1) are in line with Makino et al., (2003), who reported a considerably higher photosynthetic N use efficiency under saturated CO 2 conditions in C 4 than in C 3 leaves. Fifthly, experimental evidence suggests that Φ 2LL responds to temperature (Bernacchi et al., 2003;Yin et al., 2014). Due to the lack of understanding of this response, we empirically express the factor for the temperature effect, using a normal distribution alike equation (June et al., 2004): Finally, Equations 14, 15, and 17 require T l , and T l is solved from coupled modelling of leaf photosynthesis and transpiration: the algorithms in Supplementary Text 1 solve A and g s simultaneously; the obtained g s is used as input to the Penman-Monteith equation (Monteith, 1973;Goudriaan and van Laar, 1994) to solve leaf transpiration and T l . This procedure involves iterations, in which T l is initially set to be the same as the air temperature and then the solved T l is used for re-calculating A, g s , and leaf transpiration (Yin and van Laar, 2005).

Revising the C 4 model for simulating the cyanobacterial CCM
The single-cell C 4 photosynthesis model of von Caemmerer and Furbank (2003; see also Supplementary Text 2) can be used for simulating cyanobacterial photosynthesis . However, this model is hard to solve once it is coupled to a g s model (Equation 2 or 11). We therefore revise the above C 4 model to simulate the cyanobacterial CCM, based on the model concept of Price et al. (2011). These revisions are: (i) set g bs to a high value to mimic g ch (conductance of the chloroplast envelope to CO 2 ); (ii) set V p as if it stands for the combined rate of cyanobacterial bicarbonate transporters; (iii) set R m = R d ; and (iv) re-estimate f cyc and x, in view of the fact that extra ATP required for the cyanobacterial CCM also comes from the cyclic e − pathway (Shikanai, 2007). The ATP cost of bicarbonate transport may be lower than that of the C 4 CCM (Price et al., 2013;Furbank et al., 2015). Two single-gene transporters (BicA and SbtA) that have been well characterized in cyanobacteria are considered here, and Price et al. (2011) estimated that the two transporters require 0.25 and 0.50 ATP per transport event, respectively (so, φ in Equation 5 is 0.75). We re-estimated x as 0.2 and f cyc as 0.18 (Table 1)

Setting scenarios of improved leaf photosynthesis for simulation
Major routes in enhancing photosynthesis will be examined, except for the photorespiratory bypass. Modelling this bypass would require more complicated algorithms and parameters (von Caemmerer, 2013), which cannot be straightforwardly implemented to simulate field environments where modelling of g s is also needed. We examined the impact of improving g m (Tholen et al., 2012;Flexas et al., 2013), improving S c/o (Zhu et al., 2004;Parry et al., 2011), introducing the C 4 mechanisms into C 3 crops (von Caemmerer et al., 2012), and using cyanobacterial bicarbonate transporters and the CCM (Price et al., , 2013. Given that efforts to engineer these routes, especially the latter two, into new crops will most probably make progress step-wise, we propose the following nine routes (Table 2): (1) Improving g m , where the slope of g m25 versus leaf N (χ gm25 ; see Equation 16) is set from its default value 0.125 (Table 1) to be three times higher (i.e. 0.375).
(2) Improving S c/o , where S c/o25 is set from its C 3 default value 3022 (Table 1) to 4427, the observed S c/o25 for the nongreen alga Griffithsia monilis (Whitney et al., 2001). (4) Introducing C 4 biochemistry, where the C 4 photosynthesis model is used with C 4 kinetic constants (Table 1) while setting g bs as high as the C 3 default g m (i.e. setting the slope of g bs25 versus leaf N; χ gbs25 ; see Equation 16) to 0.125.
(6) Engineering the complete C 4 mechanism, where C 4 kinetic constants (Table 1) combined with a low χ gbs25 (0.007) is used in the C 4 photosynthesis model.
(7) Engineering cyanobacterial single-subunit bicarbonate transporters (BicA and SbtA), where the above revised C 4 model is combined with the default C 3 parameters with χ gbs25 = 0.125 and the revised values for φ, x, R m , and f cyc ( Table 2).
(8) Adding a more elaborate cyanobacterial CCM, whereby the carboxysome shell proteins are expressed in chloroplasts to enrich the CO 2 level around Rubisco similar to the level in the C 4 bundle-sheath compartment. This route assumes that the chloroplastic C 3 Rubisco can be reorganized into effective carboxysome structures and other requirements for carboxysome to function are optimized . So, the same model and parameter values as for Route 7 are used, except for χ gbs25 which is now set to a lower value of 0.007 as for C 4 bundle-sheath conductance.
(9) A complete cyanobacterial CCM installed. The complete cyanobacterial CCM will require replacement of the C 3 Rubisco with a cyanobacterial Rubisco in order to take advantage of better kinetic properties in a high-CO 2 carboxysome . Based on the expectation of engineering the cyanobacterial CCM that approaches photosynthetic capacities typical of C 4 plants (Price et al., 2013), we used χ Vcmax25 and χ Jmax25 of C 4 photosynthesis (Table 1) for this route. So, this route has the low ATP cost of the cyanobacterial CCM as well as a high enzymatic capacity per unit N to mimic the complete cyanobacterial CCM.
Simulation results of all nine routes will be compared with those of the default in which the C 3 photosynthesis model with the C 3 parameter values in Table 1 is used.

Modelling daily canopy photosynthesis and transpiration
In GECROS, instantaneous canopy photosynthesis and transpiration were calculated using the sun/shade model of de Pury and Farquhar (1997), in which the sunlit and shaded portions of the canopy each are considered as a big leaf, and the above leaf-level model is applied. Assuming an exponential profile of leaf N, total photosynthetically active N for each portion was calculated, and N-dependent photosynthetic parameters V cmax , J max , g m , g bs , and ε p (see Equation 16) were then scaled up accordingly to each portion of the canopy. Instantaneous rates were scaled up to daily total, using the Gaussian integration (Goudriaan, 1986) to account for any asymmetric diurnal courses of radiation and temperature. These approaches for spatial and temporal extensions apply to the case in the absence of water limitation.
In the presence of water limitation, the available water is partitioned between sunlit and shaded leaves according to the relative share of their potential transpiration (E p ) to obtain their actual transpiration (E a ). The diurnal course of available water is assumed to follow that of radiation. Based on the Penman-Monteith equation, the actual transpiration is transformed into the actual level of stomatal resistance to water vapour (r sw,a ) (Yin and van Laar, 2005): where r sw,p is the stomatal resistance to water vapour in the absence of water limitation [ = 1/(1.6g s ), where g s is solved from the algorithm in Supplementary Text 1]; r bh and r bw are the boundary-layer resistance to heat and to water vapour, respectively; γ is the psychrometric constant; and s is the slope of the saturated vapour pressure as a function of temperature (kPa °C −1 ). The actual r sw,a was converted into the actual g s , which can be used as input to the analytical quadratic model (see Supplementary Text 3) to estimate the instantaneous actual photosynthesis of the sunlit and shaded leaves. The Gaussian integration was again used to obtain the daily total of the actual photosynthesis. Equation 18 suggests that the impact of water deficit is mainly via stomatal conductance; any non-stomatal effect of water deficit is not modelled in GECROS, except when accounting for changes in T l under drought.

Crop simulation approaches
Simulations were conducted for three sites, Los Baños (14°6'N, 121°9'E; the Philippines), Nanjing (32°56'N, 118°59'E; China), and Shizukuishi (39°41'N, 140°57'E; Japan), representing tropical, subtropical, and temperate rice-growing conditions, respectively, using 31 year  baseline weather data and the present atmospheric [CO 2 ] of 400 μmol mol −1 . We also ran the model under the climate scenario for 2050, at which the expected [CO 2 ] is ~550 μmol mol −1 and air temperature is 2 °C higher than the baseline (Li et al., 2015). As we only examined the impact of changed leaf photosynthesis on crop productivity, we decoupled the GECROS soil module and used only the crop module for simulation to avoid any confounding effects from uncertainties in simulating soil processes. Potential production was simulated by setting the daily water supply to the crop as non-limiting. Water-limited production was simulated by setting the daily available water for evapotranspiration to no more than 50% of seasonal average daily transpiration simulated for the potential production, which was 1.97, 1.63, and 1.34 mm H 2 O d −1 for Los Baños, Nanjing, and Shizukuishi, respectively. Daily N supply was set in such a way that the accumulated N uptake by the crop followed the sigmoid curve of Yin et al. (2003) and that the total uptake at maturity reached 20 g N m −2 , equivalent to the N uptake in high-yielding rice experiments (Setter et al., 1994). Model parameters for phenology were calibrated (Table 3) so that simulated baseline crop duration was in line with that of the standard cultivar at each site (Li et al., 2015). Crop models are less accurate in predicting spikelet number and therefore harvest index than in predicting crop mass (Boote et al., 2013;Li et al., 2015). To minimize the impact of this uncertainty, we used the simulated total shoot mass (excluding dead leaves) at maturity as the proxy for crop productivity. Input parameters were set as the default values of GECROS for rice (Yin and van Laar, 2005) and those relevant to our study are given in Table 3. As C 4 enzyme kinetic parameters, especially their temperature responses, are less certain than the C 3 counterparts (Boyd et al., 2015), additional analysis was conducted for C 4 simulation (Supplementary Text 4). Similarly, because the exact ATP cost of bicarbonate transport is uncertain (Fridlyand et al., 1996;McGrath and Long, 2014), sensitivity analysis was conducted for Route 9 with regard to this cost (Supplementary Text 5). Further details about GECROS are given in Supplementary Text 6, and source codes of the full GECROS model can be obtained upon request.

Simulated leaf photosynthesis
All routes could increase A in the light-saturated region, especially Routes 6 and 9 (Fig. 1). In the light-limited region, the impact of the routes was smaller, and Routes 4, 5, and 6 in fact had a negative effect (see the inset of Fig. 1). This negative impact is associated with the two extra ATPs required for PEP regeneration in the C 4 cycle (von Caemmerer and Furbank, 1999), for which a high f cyc is required  Table 1). In Route 4 where this high ATP cost was not compensated by an effective CCM to suppress photorespiration, the negative effect was particularly high. Because these routes act differently for the light-saturated and limited regions, the curvature in the light response curve was diverse (Fig. 1) despite the same curvature factor θ (0.8; Table 1) used for Equation 12 describing the light response of PSII e − transport rate for all these curves.

Simulated canopy photosynthesis
Not surprisingly, the calculated daily canopy photosynthesis (A canopy,d ) increased with increasing LAI (Fig. 2), due to a higher interception of photosynthetically active radiation (PAR) at higher LAI. Also, the light response curve of A canopy,d became increasingly linear with increasing LAI, because at high LAI, leaves in the canopy are predominantly light limited, and within the light-limited range leaf photosynthesis increases almost linearly with light level (Fig. 1). Because the difference in leaf photosynthesis among the routes was mainly recognized in the light-saturated region (Fig. 1), the ratio of A canopy,d of photosynthesis-enhancing routes to that of the default C 3 route increased with increasing radiation level, and decreased with increasing LAI (Fig. 2). A canopy,d of Route 4, compared with the default C 3 route, was notably lower, regardless of the radiation level, when LAI was ≥3 (Fig. 2).

Default simulation for crop durations and mass production
Using GECROS, we simulated crop duration and mass production. A 2 °C warming for 2050, relative to the present climate, was simulated to shorten crop duration by ~5, 10, and 20 d, for tropical, subtropical, and temperate environments, respectively (Table 4). This different effect across the environments is due to the fact that temperature during the growing season in the tropics is around the optimum value, at which warming is expected to have a smaller effect than at the other sites where the growing season temperature is mostly in the range where development rate increases greatly with warming.
Despite the shorter duration, simulated aboveground mass at crop maturity increased for 2050 compared with the present climate (Table 4), largely due to CO 2 elevation from 400 μmol mol −1 to 550 μmol mol −1 . This is because we implicitly a One thermal day is equivalent to one calendar day if the temperature at each moment of the day is always at the optimum. b Values used for Los Baños (the Philippines), Nanjing (China), and Shizukuishi (Japan), respectively. The STTIME value for Los Baños is for the dry season there (which is the season with the high yield potential), and that for Nanjing is for single-cropping rice (that is predominant in the region, compared with the double-cropping rice where rice is planted twice per year). assumed that future breeding can develop rice cultivars capable of coping with any effect of warming on spikelet sterility, so the effect of CO 2 elevation was dominant. However, a recent FACE (free-air CO 2 enrichment) study (Cai et al., 2016) using a present cultivar showed that yields of rice were decreased by 17-35% under the combination of elevated CO 2 and temperature, compared with the ambient condition, due to fewer filled grains at the elevated temperature. As expected, water limitation decreased mass production, but increased water use efficiency (WUE) ( Table 4). The WUE differed little among the three sites, but was higher for 2050 than for the present climate, partly due to increased A canopy,d and partly due to generally decreased canopy transpiration under the 2050 climate (Table 4). The reduced canopy transpiration for the 2050 climate was largely a result of partial stomatal closure induced by higher [CO 2 ] (e.g. Wong et al., 1985).

Impact of photosynthesis-enhancing routes on crop mass production
Compared with the default C 3 photosynthesis, all routes increased aboveground mass production, except for three cases for the potential production-2050 climate combination where the benefit from Route 5 was virtually nil or slightly negative (Table 5). In general, the benefit from Routes 1, 2, and 7 was ≤10%. All routes resulted in higher benefits in the presence of drought than in the absence of drought, and under the present climate than for 2050. This could be explained by the shape of a diminishing return for A-C i curves, because drought and the present climate both result in a lower C i compared with the potential production level and the 2050 climate, respectively. Route combinations had an equal effect to, or a larger effect than, the sum of the routes acting alone. For example, the benefit from Route 3 (improving both g m and S c/o ) was about the sum of the benefits from Route 1 (improving g m ) and Route 2 (improving S c/o ) for the potential production and was higher than the sum of the two for the water-limited condition (Table 5).
The benefit from Route 6 (the complete C 4 mechanism) was considerably higher than the sum of the benefits from Route 4 (C 4 biochemistry components) and Route 5 (Kranz anatomy components for low g bs ) for any condition (Table 5). This result suggests that the ongoing programme of installing C 4 photosynthesis into C 3 crops (von Caemmerer et al., 2012), if successful, needs to engineer the complete C 4 mechanism. The benefit of a partial engineering is only marginal or even counter-productive under future high-CO 2 environments (see Route 5 in Table 5), because of the high ATP cost for operating the C 4 cycle. In an earlier preliminary simulation analysis Fig. 2. Calculated daily canopy photosynthesis of the default C 3 (0) and nine (1-9) photosynthesis-enhancing routes in response to daily incoming global solar radiation, for four different sizes of canopy [leaf area index (LAI) = 1, 3, 5, and 7, respectively]. The calculation was made using the model described in the text, based on parameter values listed in Tables 1 and 2 and the following input conditions: C a = 400 μmol mol −1 , T l = 25 °C, canopy average leaf nitrogen = 2.3 g m −2 , n b = 0.3 g m −2 , VPD = 2.0 kPa, daylength = 13 h d −1 , fraction of diffuse irradiance = 0.2, and canopy average leaf angle (from horizontal) = 65 °. Light extinction coefficient and nitrogen extinction coefficient required for the canopy photosynthesis model were calculated using the formulae as in GECROS with leaf scattering coefficient of 0.2 for PAR. (Yin and Struik, 2008), we showed that a low g bs alone would increase rice yield in the tropics by ~25%. However, that analysis used an arbitrarily low g bs and a version of the C 4 model that assumes an H + :ATP ratio of 3, whereas a recent analysis suggested that this ratio is most probably 4 , suggesting that the model version Yin and Struik (2008) used may have underestimated quantum requirement for the C 4 CCM. From our present analysis, even with the complete C 4 mechanism, its advantage over the C 3 default was simulated to be >50% only under the combination of water limitation and the present climate; for other conditions, its advantage ranged between 22% and 40% ( Table 5).
As engineering for the complete Kranz anatomy is challenging, the CCM in cyanobacteria, which is probably less expensive energetically, has been suggested as an obvious alternative to engineer (Price et al., , 2013Furbank et al., 2015). Our simulation showed that the simplest form for the cyanobacterial CCM with bicarbonate transporters, Route 7, had a marginal advantage (Table 5). Pengelly et al., (2014) showed that tobacco plants transformed with the BicA transporter had no discernible effect on CO 2 assimilation rates, suggesting that BicA was either not located or not activated correctly. Our simulation (Table 5) showed that a more elaborate cyanobacterial CCM, where the carboxysome shell proteins were expressed to enrich the CO 2 level around Rubisco in chloroplasts (Route 8), had a higher advantage than the equivalent C 4 CCM (Route 5), largely because of a lower ATP cost assumed for the cyanobacterial CCM. However, its benefit was lower than from the complete C 4 mechanism, namely Route 6, which includes the additional mechanism that C 4 plants have a considerably high carboxylation rate per unit leaf N (Evans and von Caemmerer, 2000;Makino et al., 2003). The complete cyanobacterial mechanism (Route 9), which has the low energy cost for the CCM as well as the high carboxylation and etransport capacity per unit N (presumably as high as for C 4 plants), was the only route that could bring an advantage of ≥50% under any environmental conditions (Table 5). However, as the exact ATP cost for the cyanobacterial CCM is uncertain, the simulated benefits of Routes 7-9 should be considered as tentative and their real benefits might be lower (Supplementary Text 5).

Effects of enhanced leaf photosynthesis on some other crop traits
The benefit from all routes for simulated season-long canopy photosynthesis (A canopy,s ), shown in the upper rows of Fig. 3, differed from that shown in Table 5 for aboveground mass. This difference suggests that other traits were also affected by altered photosynthesis. Aboveground mass can be calculated as: where PAR int is the season-long intercepted PAR by the crop (MJ m −2 ), PLUE is the overall photosynthetic light use efficiency (= A canopy,s /PAR int , g CO 2 MJ −1 PAR), RESP is season-long crop respiration (g CO 2 m −2 ), F root is fraction of mass for roots,  -, simulations assumed that drought had no impact on phenological development, so the predicted phenology was the same under water-limited as under the potential production level.
a Water use efficiency, defined as total crop mass production divided by the amount of water transpired during the growth season. and SENES is aboveground mass lost after leaf senescence (g DM m −2 ). We obtained the data for these five traits from the GECROS output, and calculated the percentage change of each of these traits for each route relative to the C 3 default (Fig. 3).
The enhanced leaf photosynthesis from most routes also resulted in increased PAR int , although to a small extent only, up to a maximum of 8.3% (Fig. 3). As a result, the percentage increase in PLUE by the routes was generally slightly lower than that in A canopy,s (Fig. 3). The increased PAR int stemmed from an increased LAI in the early growth phase (results not shown), in line with the recent result of  showing that increased leaf photosynthesis also resulted in increased leaf area. The increased canopy photosynthesis, on the one hand, increased the component of growth respiration; on the other hand, it decreased the component of maintenance respiration which is modelled in GECROS dependent on crop N status. The net result is that most routes decreased RESP (Fig. 3). The simulated F root generally became higher with photosynthesis-enhancing routes (Fig. 3), as expected from the classical functional equilibrium theory (Brouwer, 1983). However, in the presence of water limitation, F root could be lower compared with the default values, probably because the crop maintained a comparatively high N status under drought. Because a higher photosynthesis resulted in a lower N:C ratio in the crop and GECROS modelled leaf senescence depending on the relative magnitude of N-and C-determined LAI, the most effective enhancing routes, such as Routes 6 and 9, resulted in increased leaf senescence (Fig. 3). This simulation result is in analogy to the faster senescence and reduced LAI in later stages of development for C 3 crops grown under elevated [CO 2 ] in FACE experiments (e.g. Kim et al., 2003). Sinclair et al., (2004) simulated that a 33% increase in leaf photosynthesis may translate into only a 5% increase in soybean grain yield, or a -6% change in grain yield in the absence of additional N, presumably associated with more leaf senescence.  reported experimentally a lower than expected crop yield stimulation with rising [CO 2 ].
Next, we evaluated the extent to which these secondary effects also contributed to differences in the simulated aboveground mass. This was done from the difference in the significance level of individual terms of multiple linear regression of aboveground mass versus PAR int , PLUE, RESP, F root , and SENES (Table 6). Although the decreasing effect of SENES on mass could not be identified (because of the collinearity between mass and SENES), the significant effect of the other four terms (PLUE, PAR int , RESP, and F root ) was well estimated. While the primary PLUE had the strongest effect under both potential and water-limited conditions, the secondary PAR int always had the second strongest effect. These results on the importance of the secondary effects on croplevel traits suggest that most existing simulation studies on assessing the impact of engineering photosynthetic targets are incomplete, because computation was only done for leaf photosynthesis (e.g. McGrath and Long, 2014) or for canopy photosynthesis (Zhu et al., 2004;Song et al., 2013).

Assessing the importance of individual biochemical targets
While secondary traits were affected, after all one would assess how the primary PLUE is affected by individual biochemical targets or parameters of photosynthesis. We regressed PLUE against individual photosynthetic parameters, with site included as covariate in the regression to remove any effect of possible site differences.
The regression analysis based on simulation results using the C 3 model (Routes 1-3 plus the default) indicated that manipulating S c/o affected PLUE more than manipulating g m (results not shown), consistent with the result that the percentage change in PLUE by Route 2 was higher than that by Route 1 (Fig. 3). However, the relative impact of manipulating S c/o and g m depends on the extent to which they could actually be changed. Furthermore, an improvement in S c/o may be at the cost of decreasing V cmax , because of the often observed negative correlation between Rubisco S c/o and carboxylase turnover rate (e.g. Kubien et al., 2008;Perdomo et al., 2015). This negative correlation was not considered here, in view of the fact that the non-green alga G. monilis has a high S c/o25 while maintaining a Rubisco turnover rate similar to C 3 plants (Whitney et al., 2001). The impact of the trade-off between Rubisco S c/o and carboxylase turnover on canopy photosynthesis was analysed by Zhu et al., (2004). Our simulations using the C 4 model (Routes 4-9) involved changes in values of a set of parameters. The most important ones are χ gbs25 (which determines the effectiveness of the CCM), φ (extra ATP requirement for the CCM, which determines the required f cyc and, therefore, light-limited photosynthetic efficiency), and χ Vcmax25 or χ Jmax25 (which determine light-saturated photosynthetic capacity). Other parameters (e.g. some C 3 parameters used for Route 5) had little impact on the shape and values of light response curves (results not shown). We therefore conducted the analysis of regressing PLUE versus χ gbs25 , 3 + φ (total ATP requirement per mol CO 2 assimilated, ATP req ), and χ Jmax25 (Table 7). There was little effect of sites on PLUE. All three parameters were important for any production level-climate combination. Comparatively, the CCM parameter χ gbs25 became most important for water-limited production, because a more effective CCM to elevate the CO 2 level around Rubisco can more effectively overcome the negative effect of low C i under drought. Under the potential production, especially combined with high [CO 2 ] of the 2050 climate, the photosynthetic capacity parameter χ Jmax25 and quantum efficiency parameter ATP req were comparatively more important (Table 7). These results suggest that photosynthetic capacity, quantum efficiency, and CCM strength all need improving in order to turbocharge canopy photosynthesis. Based on natural variation of leaf photosynthesis, Gu et al., (2014a) showed that quantum efficiency parameters had even higher effects than capacity parameters on rice productivity. a Site was included as the covariate in regression, with Los Baños as the reference having a coefficient of zero. b Total ATP requirement per CO 2 assimilated (= 3 + φ), i.e. 5 for C 4 photosynthesis and 3.75 for cyanobacterial photosynthesis (see the text). The five component traits are: PAR int , season-long intercepted PAR; PLUE, overall photosynthetic light use efficiency defined as season-long canopy photosynthesis divided by PAR int ; RESP, season-long crop respiration; F root , fraction of mass for roots; SENES, aboveground mass lost due to leaf senescence. Linear regression is given as:

Concluding remarks
We simulated the likely impact of major routes in ongoing programmes using transgenic technology to improve photosynthesis (Table 2). Our analysis showed that improving leaf photosynthesis can result in an increased rice mass production to a different extent (Table 5), thereby also resulting in different improvements in resource use efficiency such as WUE (Fig. 3). However, to supercharge photosynthesis significantly, engineering for a single improvement route can hardly be effective. Some single routes may be counter-productive at the canopy level. For example, installing C 4 biochemistry, if not combined with an effective CCM, is only beneficial for upper leaves of the canopy, while it has no or even a negative impact for lower shaded leaves because such a mechanism requires extra ATP for the C 4 cycle. Note that the standard C 4 model of von Caemmerer and Furbank (1999) for e − transport limitation does not explicitly consider the increased cyclic e − transport due to the extra ATP costs relative to C 3 photosynthesis, and, therefore, cannot recognize the little advantage of C 4 photosynthesis under shade. Similarly, the simulation by McGrath and Long (2014) in assessing the potential of cyanobacterial CCM took no account of the extra ATP required by bicarbonate transporters. Our simulation also showed that manipulating photosynthesis may result in unwanted secondary effects on some traits at crop level (e.g. inducing faster senescence if nutrient uptake is not increased). Therefore, the beneficial effect of the single route for high photosynthesis on increasing crop productivity may have previously been overestimated. To supercharge crop productivity, combined routes for improved CCM, photosynthetic capacity, and quantum efficiency are required.

Supplementary Data
Supplementary data are available at JXB online. Fig. S1. Simulated versus observed relationships between A and g s . Text 1. Analytical solution to the cubic equation as a result of combined stomatal conductance, CO 2 diffusion, and biochemical leaf photosynthesis models.
Text 2. Revising the C 4 photosynthesis model for simulating the cyanobacterial CCM.
Text 3. Analytical solution to the quadratic equation as a result of combined CO 2 Text 4. Analysis with respect to temperature response parameters of C 4 enzyme kinetics.
Text 5. Sensitivity analysis with respect to ATP cost for the cyanobacterial CCM.