Rubisco substitutions predicted to enhance crop performance through carbon uptake modelling

A state-of-the-art canopy modelling framework was used to predict changes in crop carbon uptake after Rubisco substitutions. We provide a new perspective and improve predictions from previous Rubisco modelling studies at the leaf level.


Canopy temperature
Leaf temperature is solved using surface fluxes of sensible heat flux , latent heat flux and net radiation at each time point (Bonan, 2019).
Where and are subtracted from the total net radiation of the site surface. The remaining energy retained by the surface is represented using a leaf heat capacity . is a function of leaf mass per unit area of a leaf. The fraction of which is water is assumed (e.g. 80% in crop species) and the fraction of which is dry matter is assumed to be 1 − . Heat capacity constants of dry matter and water are also assumed for the dry and water fractions.
is represented by the following equation: Liou (2002) provided an analytical solution for the two-stream approximation as described by Bonan (2019), (Bonan et al., 2011).
Equation (3) gives the upward flux of direct beam solar radiation above the leaves and equation (4) gives the downward flux below the leaves at depth equal to the cumulative leaf area index; ∆ =LAI. In this study, clumping index was ignored (Ω=1). Inputs for equations (3) and (4) , ↓ is the total downward direct beam solar radiation flux which can be obtained from eddy covariance data sets for each site. ωℓ = Ρℓ + τℓ is a scattering beam coefficient obtained by adding a known leaf reflectance Pℓ and transmittance τℓ.
is a known soil albedo. The proportion of direct (β 0 and diffuse beam β radiation scattered is then obtained using the following equations as described by Sellers (1985): Θℓ is the mean leaf inclination angle. cos 2Θℓ was approximated using the Ross-index ℓ which describes the departure of leaf angle from a spherical distribution. For crops ℓ was obtained from the CLM4.5 as follows: The coefficient which is the direct beam extinction coefficient (defined as the optical depth per unit leaf area) is probably one of the most important parameters in canopy modelling as it determines the sunlit and shaded fraction of the canopy models.
is defined by the following equation: cos ℓ or ( ) or ℎ is the amount of leaf area projected in the direction of the solar beam. An approximate for cos ℓ or ℎ is obtained using the Goudriaan (1988) function ( ) and the corresponding ℓ as follows: The coefficients for equation 16 are: cos adjusts cos ℓ or ( ) or ℎ onto a horizontal plane. cos is obtained from the solar zenith angle which can easily be obtained from remote sensing data or established relationships based on location (latitude and longitude) and time (i.e. time zone month, year, hour day). In this study, for each half hourly observation was approximated using the 'fishmethods' R package which adopts analytical equations by Frouin et al. (1989).
Because diffuse radiation is assumed to come from all parts of the sky (unlike direct beam radiation), the extinction coefficient for direct beam radiation was approximated by dividing the sky hemisphere into three 30 o increments and summing the transmittance obtained from these three angles to give an approximate of diffuse transmittance as follows: is the coefficient for each 30 o increment obtained using equation (15). Using the trigonometric relationship 2 sin cos ∆ gives the total contribution of each angle to total sky irradiance. Each 30 o increment of the sky hemisphere was assumed to be equally spaced apart (i.e. ∆ =0.1745 radians) (Bonan, 2019). The final coefficient for diffuse radiation is given using the following equation: Finally, the single scattering albedo ( ) for equation (13) was obtained using the following equation: Where ( ), ∅2 and ∅1 are obtained from the Goudriaan (1988) function (16).
The overall radiative balance for direct beam radiation is given by considering the upward downward flux (per leaf area) at a depth equal to the cumulative LAI : Similar equations pertain for diffuse beam radiation → ( ) but involve dropping direct beam terms in equation (3) and (4) and changing coefficients 1 and 2 (see pg.252 Bonan (2019)).
The sum of the diffuse → ( ) and scattered beam radiation → ( ) gives the radiation absorbed by the shaded leaves → ℎ (per shaded leaf area) at the cumulative leaf area index : Where the → ( ) is: The shaded radiation absorbed by sunlit leaves (per sunlit leaf area) is: Integrated over the sunlit and shaded fraction of the canopy ℎ : The sunlit and shaded LAI fractions are given as follows:

Photosynthesis model and CO2 diffusion
The Farquhar et al. (1980) model for C3 photosynthesis and Collatz et al. (1992) model for C4 photosynthesis were implemented in this study. Both models depend on a set of parameters (Table S1). Briefly, leaf net photosynthesis is the minimum of RuBP carboxylate (Rubisco) limited rate , light limited rate and the phosphate limited rate for C3 photosynthesis or phosphoenolpyruvate (PEP) carboxylase limited rate for C4 photosynthesis minus dark respiration (29).
= min( , , ) − The RuBP carboxylate limited rate is given using the hyperbolic function for C3 photosynthesis and for C4: The light limited rate is given using another hyperbolic function for C3 photosynthesis and in relation to quantum efficiency for C4: The phosphate limited rate is given for C3 and PEP carboxylase limited rate for C4 plants as follows: In equations 30-31, , , / , and are the Michaelis-Menten kinetics for the maximum velocity of Rubisco and its corresponding light limited rate , specificity for carbon dioxide (CO2) to oxygen (O2) / , Michaelis-Menten constant for CO2 , and Michaelis-Menten constant for O2 , respectively.
is the amount of atmospheric CO2 that is assumed to be in the leaf.
is the amount of atmospheric oxygen assumed to be present in the leaf. For simplicity, this is kept constant at = 210 2 −1 . Γ * is the CO2 compensation point given as Γ * = 0.5( / ). The quantum efficiency for C4 plants is assumed to be = 0.05 2 −1 ℎ . is the absorbed PAR which is obtained from equation 27 and 28 for sunlit and shaded fractions of the canopy. The triose phosphate utilization rate is obtained using the relationship: = 0.167 . Similarly, the initial slope of the CO2 response curve for PEP carboxylase limited rate for C4 plants is given: For C3 plants is obtained using the relationship: = 1.65 . Jmax is also dependent on PAR so requires adjustment for variations in PAR. In the CLM, this is given using the quadratic function and the solution is the smaller root: Θ is a curvature parameter which is set as 0.7 in this study. Ι is the amount of light used in photosystem II given as follows: Φ is the quantum yield of photosystem II which is assumed to be Φ = 0.85 and 0.5 partitions photons between both photosystems.
In this study, the co-limitation as described by Collatz et al. (1992) is used which smooths the transition between , using the following quadratic functions: The first quadratic function in equation (35) obtains an intermediate gross CO2 assimilation (i.e. between and rates) and then and are used to obtain the final gross CO2 assimilation. Θ and Θ are smoothing parameters for C3 or C4 plants.

Temperature adjustments for Rubisco kinetics
The Rubisco kinetic parameters , , Γ * , and all depend on temperature. The kinetics measured at 25 O C were adjusted using the Arrhenius equation (Bonan, 2019, Bonan et al., 2011 (36) is the species-specific heat activation value (KJ mol -1 ), is the universal gas constant, is the canopy temperature obtained from equation (1) and 25 is a species-specific kinetic parameter at 25 O C.
is the temperature adjusted kinetic parameter. Since some parameters (e.g. , , and ) vary with for simplicity some parameters were adjusted according to temperature changes of except Jmax. It is understood that the temperature response of Jmax is closely related to the lipid composition and thylakoid composition (Von Caemmerer, 2000). For simplicity, an average C3 Ha value was used for all crop species Jmax.
The temperature adjusted parameters were additionally adjusted for enzyme breakdown at extreme temperatures. Separate temperature functions with constants for C3 and C4 crops are used as it is unknown whether there are significant inter-species differences.
For C3 photosynthesis the temperature deactivation function is: Where ∆ (KJ mol -1 ) is an entropy term and is the heat deactivation value (KJ mol -1 ) ( Table  S1).
For C4 photosynthesis the temperature deactivation function is (Sellers et al., 1992: Where the parameters 1, 2, 3, 4 (K -1 ) adjust the rate of increase/decrease of at low or high temperatures.

Coupling photosynthesis models to CO2 diffusion model
Both photosynthesis models must be coupled to a model which quantifies the amount of that reaches the site of carboxylation in the leaf. Mathematical expression for full CO2 diffusion from the atmosphere to the site of carboxylation which is in the chloroplast is as follows: A is the gross photosynthesis assimilation and , , and are the CO2 conductance of the leaf boundary layer, stomata and mesophyll, respectively given as mol H2O m -2 s -1 . In Earth system models (ESMs) CO2 diffusion into the mesophyll is often ignored as it remains uncertain how to correctly parameterise in ESMs (Rogers et al., 2017). is dropped from equation (27) and the intercellular CO2 is assumed to be a product of and as follows: In equation (39) and are calculated for H2O and 1.4 and 1.6 in equation (40) adjusts the conductance for H2O to that of CO2 which accounts for the lower diffusivity of CO2 compared to H2O.
In this study was calculated using the semi-empirical Ball & Berry stomatal model as follows (Bonan, 2019, Franks et al., 2017: The Ball & Berry model assumes that stomatal conductance varies linearly with , relative humidity ℎ and surface CO2 . 0 and 1 are known minimum and maximum values. is a soil moisture stress parameter which adjusts the relationship for poorly watered conditions. This is required since 0 and 1 measurements are often done in well-watered conditions. For simplicity it was assumed that leaf ℎ was the same as surface humidity reported by the flux data sets. is obtained from using equations derived from engineering studies as follows (see p.159 Bonan (2019)): (42) and are parameters that describe the effect of wind speed on a leaf's length ℓ. In this study, ℓ was assumed to be 0.05 meters.
in equation (41) was calculated using the reported soil water content for each half hourly observation for each site as follows: is the field wilting point which is assumed to be 0.1 (or 10%). is when soil moisture is assumed to be saturated. In this study, values for different soil textures were derived from Campbell (1974). is the soil water content from each half hourly observation.
Because equation (40) is also dependent on photosynthesis and the solution to equation (40) is also an input into the photosynthesis models. Equation (40) requires a numerical approach. In this study, the Newton-Raphson root finding method was used as described by Sun et al. (2012). The Newton-Raphson approach is as follows: The initial guess for Ci is obtained using a ratio of to conversion for C3 and C4 plants. The final iteration function is: Newton Raphson method requires the derivative ′( ) of ( ). Sun et al. (2012) used a simple numerical approximate using a finite difference: This initial guess is passed through a series of equations: = min ( ( ), ( ), ) The full iteration scheme was as follows: 1. An initial guess for was obtained (47). 2. was obtained using equation (48). 3.
was obtained using a quadratic function of the Ball & Barry equation (50). 5. Final was obtained using equation (51). 6. The iterative function ( ) = ℎ( ) − ( ) (45) was formed using the initial guess from step 1 and final guess from step 5. 7. A finite difference for the derivative ′( ) was obtained by adding the initial and final . 8. The finite difference was then used to obtain an intermediate , , and final (48-51).
9. An approximate of the derivative was obtained according to equation (46) as follows; An intermediate iteration function (45) (46) and iteration function from step 6 formed the denominator.
10. Steps 2-9 were repeated until the convergence criteria was met (<0.001 µmol CO2) and the upper limit of iterations was set to N= 20.
If the upper limit was reached because convergence was not met, the bisection method was used which guarantees convergence (Burden, 2016). The bisection method continuously bisects an interval containing the correct solution 'x' until convergence is achieved.
For each iteration, the bisection method is in the form of: is the midpoint of the interval and . The full iteration scheme was as follows: 1. To obtain the first guess, , the lower interval was assumed as the lowest possible that could be achieved; Γ * . was assumed as the highest that could theoretically be achieved: (52). 2. An initial was obtained using the initial guess (48). 3. The initial was used to obtain an initial (49). 4. The initial was used to obtain an initial using a quadratic function of the Ball and Berry equation (50). 5. The initial was used to obtain an initial (51). 6. If the difference − was below the convergence criteria (<0.001 µmol CO2) then the iteration stopped, and convergence was met. 7. If the convergence criteria were not met, bracketing of the interval continued as follows: the initial ′ ′ was re-calculated using the function (48-51) to obtain ( ) and then the difference ( ) − . If the sign of the difference − and the difference ( ) − were positive, the initial midpoint c would become the next iteration otherwise if either of the signs were negative the initial remained unchanged. Similarly, if the sign of the difference − and the difference ( ) − were positive, the initial c midpoint would become the next iteration ′ ′, or remain unchanged. 8. Finally, if the signs of the initial difference − and the difference ( ) − were positive the next iteration's difference ( ) − would become − otherwise remain unchanged. Steps 2 to 8 were repeated until the convergence criteria was met.

Gradient of Rubisco kinetics throughout the canopy
Rubisco kinetics including the / , and were assumed to be the same throughout the canopy. An approximate of 25 measured at 25 O C for the top of the canopy was obtained as described by Houborg et al. (2013): 25 is the maximum catalytic rate of Rubisco measured at 25 O C, is the fraction of leaf nitrogen in Rubisco, leaf nitrogen content and 0 is a conversion factor. was obtained as follows: : is a carbon to nitrogen ratio. is the specific leaf area at the top of the canopy.
25 for the sunlit and shaded fractions of the canopy was obtained as follows (Bonan, 2019, Bonan et al., 2011:  (19) and the nitrogen coefficient that describes the vertical gradient of nitrogen was assumed to be the same for all crops ( =0.3).

Modelling seasonal leaf area index
LAI is a major driving force of the simulations (Bonan, 1993). The CLM uses a complex carbonnitrogen biogeochemical model to predict seasonal LAI prospectively. Since the aim in this study was to retrospectively model changes in photosynthesis, there are less complicated methods for obtaining seasonal trends of LAI for retrospective analysis.
To rule out any errors in the simulations arising from inaccurate LAI estimates, three LAI retrieval methods were tested for their ability to accurately simulate flux site net CO2 assimilation trends retrospectively. The three LAI methods included a normalised vegetation index (NDVI) method derived from flux tower incoming and outgoing photosynthetically active radiation (see Wilson and Meyers (2007)), eight day MODIS satellite LAI product (https://modis.ornl.gov/sites/) and the growing production day (GPD) method which uses estimates of gross primary productivity (GPP) as a proxy of canopy size but modified here to use tower derived GPP rather than estimates from the MODIS GPP algorithm (see Xin et al. (2019)). 2005 growing season of the US-Bondville maize site was chosen for the simulations.  (Sharwood et al., 2016, Hermida-Carrera et al., 2016.