The role of chloroplast movement in C4 photosynthesis: a theoretical analysis using a three-dimensional reaction–diffusion model for maize

A novel 3D reaction–diffusion model of gas transport during C4photosynthesis revealed that the movement of mesophyll chloroplasts increased rate of photosynthesis at the cost of increased leakiness.


Introduction
Photosynthesis in C 4 plants is empowered by the CO 2 concentrating mechanism (CCM) that enables a higher CO 2 concentration in bundle sheath chloroplasts, allowing ribulose-1,5-bisphosphate carboxylase/oxygenase (Rubisco) to operate near a CO 2 saturation level. During C 4 photosynthesis, the CO 2 that diffuses from the intercellular air spaces (IAS) to the mesophyll cytosol is hydrated to bicarbonate, which reacts with phosphoenolpyruvate (PEP) in the presence of phosphoenolpyruvate carboxylase (PEPC) to form C 4 acids. The C 4 acids diffuse to bundle sheath cells and are decarboxylated to release CO 2 for fixation by Rubisco in bundle sheath chloroplasts (Furbank and Hatch, 1987;Furbank, 2011). The affinity of PEPC for bicarbonate is much higher than that of Rubisco for CO 2 , and the PEP carboxylation is much faster than RuBP carboxylation by Rubisco. This increases the CO 2 concentration in bundle sheath chloroplasts, which supercharges the catalytic sites of Rubisco, reduces photorespiration, and thus enhances the photosynthetic rate (Kanai and Edwards, 1999). An efficient operation of the CCM requires the coordinated functioning of mesophyll and bundle sheath cells (Pengelly et al., 2010;Kromdijk et al., 2010Kromdijk et al., , 2014Sun et al., 2012Sun et al., , 2014Bellasio and Lundgren, 2016) as the CCM is more costly in ATP to regenerate PEP in mesophyll cells. The concentric arrangement of mesophyll, bundle sheath, and vasculature, the relative size of mesophyll and bundle sheath cells, and chloroplast abundance and distribution impact the light harvesting potential of bundle sheath cells Stata et al., 2014;Bellasio and Lundgren, 2016) influencing the operation of the CCM.
Mesophyll chloroplasts of C 4 plants are randomly distributed along the cell walls while those of bundle sheath cells are typically close to mesophyll cells (centrifugal location), which is an arrangement common in the classical NADP-ME subtype of C 4 plants that includes major C 4 crops (Dengler and Nelson, 1999). Upon exposure to high light for several hours or at extremely high light intensity at midday, mesophyll chloroplasts of finger millet (Eleusine coracana, a NAD-ME species) exhibited movement towards bundle sheath cells (referred to as aggregative movement) (Maai et al., 2020a). Under similar conditions, mesophyll chloroplast of NADP-ME species, like sorghum and maize, aligned between veins, along the sides of anticlinal walls, and parallel to the direction of light to allow light to pass through (referred to as avoidance movement) (Yamada et al., 2009;Maai et al., 2020b). When exposed to blue light after dark adaptation for several hours, mesophyll chloroplasts of maize and sorghum showed aggregative movement (Ryu et al., 2014). In addition, the aggregative movement has recently been shown to occur in a number of other C 4 plants (Kato et al., 2022). Unlike mesophyll chloroplasts, the intracellular arrangement of bundle sheath chloroplasts was unaffected in finger millet, maize, or sorghum (Taniguchi et al., 2003;Yamada et al., 2009;Maai et al., 2020a). C 4 crops experience high light for a prolonged time in their natural environment, especially at the top of the canopy, which may compromise yield (Anderson et al., 2021;Wang et al., 2022). Given the variety of arrangements of mesophyll chloroplasts in response to high light intensity, it is important to understand what the consequences of each chloroplast configuration are for the leaf energetics and net photosynthetic rate (A n ). The analysis requires knowing the amount of light absorbed by mesophyll and bundle sheath chloroplasts to analyse the cell-specific ATP and NADPH demand and production, which are hard to measure but important parameters of C 4 energetics (Bellasio and Lundgren, 2016;Yin and Struik, 2018). The CCM efficiency is quantified by leakiness (Φ), defined as the ratio of CO 2 leakage flux from bundle sheath cells to the rate of PEP carboxylation (Kanai and Edwards, 1999;Yin and Struik, 2018).
Chloroplast movement in C 3 plants in response to light and its effect on CO 2 assimilation and productivity have been extensively investigated (Tholen et al., 2008;Wada, 2013;Gotoh et al., 2018). In C 4 plants, chloroplast movement following short term changes in light quality in maize was proposed to result in better coordination between CO 2 fixation in mesophyll and bundle sheath cells, and thus reduced Φ (Sun et al., 2014). The aggregative movement is assumed to boost photosynthesis through improved light penetration or by shortening the diffusion path of metabolites and re-fixation of leakage CO 2 (Taniguchi et al., 2003;Yamada et al., 2009;Maai et al., 2011Maai et al., , 2020aTaniguchi and Cousins, 2018). However, it was found that the extent of aggregative movement did not differ for different bundle sheath chloroplast positions and for the presence or absence of a suberin lamella in bundle sheath cells, factors that may affect CO 2 leakage, its refixation and metabolite transport (Kato et al., 2022). By contrast, the avoidance movement, which is associated with photoprotection of mesophyll chloroplasts in sorghum, may decrease A n (Maai et al., 2020b).
The effect of chloroplast movement on the energetics of C 4 plants and photosynthetic rate can be analysed in a controlled manner using modeling due to otherwise simultaneous changes in leaf biochemistry, stomatal conductance, and anatomy of C 4 plants upon exposure to high light stress (Maai et al., 2020b;Anderson et al., 2021). C 4 leaf energetics were recently examined using biochemical models and flux balance analysis (Bellasio and Lundgren, 2016;Bellasio, 2017;Struik, 2018, 2021;Bellasio and Farquhar, 2019). Some of these studies (Bellasio and Lundgren, 2016;Yin and Struik, 2018) accounted for the effect of leaf anatomy on light propagation indirectly using the Beer-Lambert law (Bellasio and Lundgren, 2016;Yin and Struik, 2018). The light propagation was modeled using the meshed Monte Carlo method, which, unlike the classical Monte Carlo method for photon transport, uses tetrahedral meshes for better representation of the structural intricacies of biological tissues (Watté et al., 2015;Ho et al., 2016). Reaction-diffusion (R-D) models of CO 2 transport account for the effect of leaf anatomy on light and CO 2 distribution directly . R-D models were applied to analyse gas exchange during photosynthesis using two-dimensional (2D) (Ho et al., 2012;Retta et al., 2016a;Berghuijs et al., 2019) or three-dimensional (3D) (Tholen and Zhu, 2011;Ho et al., 2016;Wang et al., 2017;Xiao et al., 2023) leaf anatomy. However, 3D R-D models provide the capacity for more accurate and mechanistic analysis of CO 2 and light propagation as emphasized recently (Earles et al., 2019). Exploring the influence of a 3D arrangement of mesophyll and bundle sheath chloroplasts on photosynthesis requires a realistic 3D geometry at tissue level for accurate modeling of light propagation and diffusion of gases for the irregular geometry of the IASs Xiao et al., 2016;Théroux-Rancourt et al., 2017;Earles et al., 2019). Yet the application of such models to C 4 photosynthesis has received less attention than for C 3 plants.
To quantify the effect of chloroplast movement on C 4 energetics, photosynthesis, and CCM efficiency, we develop a 3D model of C 4 photosynthesis consisting of a R-D model of gas exchange coupled with an ATP and NADPH production model (Yin and Struik, 2018) and a light propagation model (Watté et al., 2015) using a 3D representation of a real maize (Zea mays L.) leaf as a model system. We hypothesized that the aggregative movement of mesophyll chloroplasts improves photosynthesis through improved light propagation, increased energy production, and reduced leakiness. The avoidance movement of mesophyll chloroplasts is hypothesized to decrease photosynthetic rate at high light intensity.

Plant sample, experiments, and reconstruction of the 3D leaf anatomy
The plant growth conditions, gas exchange, and chlorophyll fluorescence measurements were described in detail previously (Retta et al., 2016a). In brief, maize plants (hybrid 2-02R10074) were grown receiving contrasting nitrogen treatments. The growth conditions were ambient [CO 2 ] of 380 µmol mol −1 , relative humidity of 60-80%, temperature of 23 °C during the day and 16 °C during the night, and a photoperiod of 16 h. The response of net photosynthetic rate to CO 2 and light were measured in four replicates. The light response curve was measured at an ambient CO 2 level of 250 μmol mol −1 to promote photorespiratory conditions. We used these data for comparisons with the predicted responses of photosynthesis from the 3D C 4 photosynthesis model developed here.
To reconstruct the 3D geometry (referred to as the 'default' geometry), plant tissues from the same set of maize plants measured for gas exchange and chlorophyll fluorescence as cited above were imaged using light microscopy (Olympus BX-51, at ×40 magnification). Sample preparation, fixation, and imaging were reported previously (Retta et al., 2016b). A total of 120 leaf sections at 1-µm intervals were cut using one leaf from a sample that was least affected by the sample processing ( Supplementary Fig. S1A, B). The images were first aligned using the 3D image-processing software Avizo Fire (VSG, France). The images included three veins each surrounded by a layer of bundle sheath cells, each of which was surrounded by a single layer of mesophyll cells forming a repetitive concentric arrangement of bundle sheath and mesophyll cells ( Supplementary Fig. S1A). Therefore, using the symmetry of the geometry, only one concentric arrangement of mesophyll and bundle sheath cells along with the epidermis and IAS was considered by cropping the images to a representative region of interest, resulting in 3D images having a pixel resolution of 0.67 µm.
The regions of mesophyll and bundle sheath cells, chloroplasts, epidermis, IAS, and vascular bundle ( Supplementary Fig. S1A) were labelled semi-automatically. Deformations of the epidermis on some of the images (Supplementary Fig. S1A) and missing sections were corrected for by interpolating through the consecutive images in Avizo. The chloroplasts of mesophyll cells were identified clearly thanks to their random arrangement, while numerous clustered chloroplasts of bundle sheath cells were segmented as a layer of centrifugally (towards mesophyll cells) arranged chloroplasts. Other microstructural features such as mitochondria, cytosol, and vacuole could not be identified due to insufficient resolution and contrast. Thus, mitochondria were assumed to be uniformly distributed in the cytosol, and vacuoles were created in MATLAB (version 7, The MathWorks, Inc., Natick, MA, USA) assuming the vacuole occupies 60% of the cell volume and accounting for sufficient cytosol volume required for CO 2 hydration (Taiz, 1992). Two stomata identified on the light microscopy images were digitized as cylindrical pores. Two additional stomatal pores were created assuming a stomatal ratio of 1.0 so that the stomatal index matches that reported for maize (Driscoll et al., 2006;Zheng et al., 2013). The anatomical properties measured from the 3D geometry were compared with those reported specifically for maize or broadly for NADP-ME C 4 subtype plants.

Light propagation modeling and measurement of leaf optical properties
The propagation of light within the 3D geometry of maize leaf was modeled using a meshed Monte Carlo method (Watté et al., 2015). Optical properties of mesophyll and bundle sheath cells were computed based on the Mie scattering theory from the distribution of organelles in these cells (Supplementary Tables S2, S3) (Aernouts et al., 2014). The mean radius of organelles such as mitochondria, peroxisomes, nucleus, Golgi stack, and ribosome-like structures was calculated from the reported values of radius from the literature (Supplementary Tables S2). Using the mean radius and assuming spherical objects, a particle size distribution was calculated to estimate the number of organelles per unit volume. The latter was used to calculate the bulk scattering coefficients and anisotropy factor of cell media that were contributed by the optical properties of cell components such as mitochondria, peroxisomes, nucleus, Golgi stack, and ribosome-like structures (Supplementary Table S3). The scattering coefficient of chloroplasts was calculated from that of the grana. The thickness of the cell wall in our model was 160-190 nm, which is much less than the wavelength of light at which the Mie scattering theory is valid. Thus, the cell wall was not explicitly modeled but its scattering effect was accounted for by lumping the scattering coefficient of cell wall with that of the chloroplasts.
Absorption coefficents of the mesophyll and bundle sheath chloroplasts were calculated from the chlorophyll content of the leaf (Supplementary Table S3). The leaf chlorophyll content was assumed to be 539 µmol m −2 , which is an average values of 501 µmol m −2 for maize having high nitrogen content and 579 µmol m −2 , the maximum value reported for the NADP-ME subtype, for sorghum (Ghannoum et al., 2005). The content of chlorophyll in bundle sheath cells can varry from 33% to 75% of the leaf chlorophyll content (Ghannoum et al., 2005;Bellasio and Lundgren, 2016). It was assumed that the chlorophyll distribution was 50% in mesophyll and 50% in bundle sheath chloroplasts. The absorption coefficient of the epidermis, vacuole, and cytosol was set at 100 m −1 .
The propagation of blue light (470 nm) and red light (665 nm) was modeled separately by assuming the illumination to be distributed  (Frank et al., 1996;Lide, 1999;Juurola et al., 2005 Pocker and Miksch (1978), Wang et al. (2014) uniformly over the top surface of the leaf. The resulting absorption profiles for red and blue light were combined with weighting fractions of 10% blue and 90% red (corresponding to the light source used during the gas exchange measurements) to derive an overall light absorption profile (Supplementary Table S3). The resulting light profile was then used to compute the potential rate of local electron transport and ATP and NADPH production in the mesophyll and bundle sheath chloroplasts (see 'ATP and NADPH production model' below).
To validate the light propagation model, light reflectance and transmittance were measured with the integrating spheres set-up elaborated by Aernouts et al. (2014), which was used previously by Ho et al. (2016). For light reflectance and transmission measurements, a different set of maize plants from another cultivar (cultivar P8057) were used as the measurement was not carried out at the same time as the gas exchange measurement. The plants were grown in a greenhouse at the facilities of Wageningen University & Research, The Netherlands. The set point for relative humidity during growth was 70% and that for temperature was 27 °C during the day and 22 °C during the night. The photoperiod was 12 h. These conditions were like the ones used to grow the maize plants for which the earlier-mentioned gas exchange data were acquired. In addition, cultivar P8057 received a total of N, P, K, and Mg of 1020, 710, 1630, and more than 160 mg per pot, respectively. Hybrid 2-02R10074A received 500 mg-1250 mg N per pot while total P, K, and Mg were 330, 1250, and 500 mg per pot, respectively (Retta et al., 2016a). Therefore, the nitrogen applications, which influence chlorophyll content, were similar. The leaf thickness of a few randomly selected leaves was 200 ± 50 μm. Similarly, leaf thickness of hybrid 2-02R10074A ranged from 208 to 237 μm for the range of nitrogen treatments the cultivar received (Retta et al., 2016b).
Reflectance and transmittance were measured from four plants using 2-mm samples from the mid-position of the maize leaves. The samples were placed between two integrating spheres and illuminated with a super continuum laser-based tunable light source to measure the total reflectance and transmittance in the wavelength range from 550 nm to 850 nm at 5 nm intervals. The measured total reflectance and transmittance were then compared with those computed with the light propagation model.

Microscale gas exchange model
The microscale model for gas exchange of a C 4 plant leaf was developed in analogy to the 3D R-D model for C 3 leaves presented by Ho et al. (2016), and by extending a 2D R-D model of gas exchange of a C 4 leaf (Retta et al., 2016a). From the biochemical model of C 4 photosynthesis (von Caemmerer and Furbank, 1999;), we developed the equations to account for 3D fluxes of CO 2 , bicarbonate, and O 2 within the maize leaf and the rates of the associated processes. CO 2 diffusion through the stomata was modeled by adjusting the effective diffusivity of CO 2 through the stomatal pores to match the measured stomatal conductance without artificially changing stomatal geometry ( Supplementary Fig. S2). Due to uncertainties of combining detailed metabolite pool concentrations with inhibition of enzyme activity by end-products, the inter-cell-type transport was assumed to be in a steady state and non-limiting in our model, as were the gas exchange measurements. Parameters of the 3D model are given in Table 1. Symbols, definitions, and units of the physical properties and photosynthetic parameters are given in Supplementary Table S1.

Biochemical model of C 4 photosynthesis and hydration rates
The gross volumetric rate of CO 2 fixation in bundle sheath chloroplasts, V * c , is given by the minimum of Rubisco-limited carboxylation rate (first term in brackets in Eq. 1), and electron-transport-limited carboxylation rate (second term in brackets in Eq. 1) based on ATP supply according to the biochemical model of C 4 photosynthesis (von Caemmerer and Furbank, 1999;Yin and Struik, 2012): where K m,C and K m,O are Michaelis-Menten constants of Rubisco for CO 2 and O 2 , respectively; γ * = 0.5/S C/O , where S C/O is the relative CO 2 / O 2 specificity factor of Rubisco; V * c,max is the volumetric rate of maximum capacity of Rubisco carboxylation (calculated using Eq. 13); and j * ATP,C3 is the volumetric potential rate of ATP used for C 3 cycle calculated by assuming 60% of the total ATP production ( Eq. 24) is allocated to the C 3 cycle (von Caemmerer and Furbank, 1999). [CO 2 ] and [O 2 ] are the concentrations of CO 2 and O 2 in the bundle sheath chloroplasts, respectively.
The volumetric rate of PEP carboxylation (V * p ) is given as a minimum of the enzyme-limited rate (first term in brackets in Eq. 2) and the electron-transport-limited rate of PEP carboxylation (second term in brackets in Eq. 2) (von Caemmerer and Furbank, 1999;Yin and Struik, 2012;Boyd et al., 2015):  (2021), the fraction of bundle sheath surface area occupied by plasmodesmata (Sowiński et al., 2008) Relative viscosity η Cell media 2 Unitless Assumed where K h is the Michaelis-Menten constant of PEPC for bicarbonate and V * p,max is the maximum catalytic rate of PEPC (calculated using Eq. 13). The volumetric rate of ATP used for the C 4 cycle in mesophyll cells, j * ATP,C4 , was calculated by assuming 40% of the total ATP production (Eq. 24) is allocated to the C 4 cycle (von Caemmerer and Furbank, 1999).
The net hydration of CO 2 in the presence of CA, B CA, is given by (Spalding and Portis, 1985;Tholen and Zhu, 2011): where V * CA,max is the maximum catalytic activity of CA per unit leaf volume;  is the concentration of bicarbonate ions; K CO2 and K HCO3 − are the Michaelis-Menten constants of CA hydration and dehydration, respectively; K eq is the equilibrium constant for CA; and [H + ] is the concentration of H + ions.
Since bundle sheath cells of maize contain little CA, it is assumed that hydration of CO 2 proceeds non-enzymatically in these cells (Burnell & Hatch, 1988). The net non-enzymatic hydration rate of CO 2 per unit leaf volume, B NCA , is given by : where k 1 is the CO 2 hydration rate constant; k 2 is the CO 2 dehydration rate constant; K is the acid dissociation constant for H 2 CO 3 ; and [H + ] is calculated from the pH of the cell compartments such as chloroplast stroma and cytosol. Equation 4 is assumed valid everywhere in the liquid phase.
CO 2 transport inside a C 4 leaf The R-D equation for CO 2 transport at steady-state is: where ∇ is the gradient operator; D CO2 is the diffusivity of CO 2 in either liquid (cells) or IAS; G stands for source or sink terms in each cell component; R * is the volumetric rate of CO 2 release through respiration in the cytosol of mesophyll and bundle sheath cells, epidermis and vascular bundles (Eqs 8,9); A * n is the volumetric rate of net photosynthesis which includes CO 2 fixation in bundle sheath chloroplasts (Eq. 1) and CO 2 release by respiration (Eq. 9) and photorespiration in bundle sheath cytosol (Eq. 10); and V * p is the average (indicated by the over score) rate of CO 2 production through decarboxylation of C 4 acids in the bundle sheath chloroplasts per unit leaf volume (Eq. 11). As in the standard C 4 photosynthesis model (von Caemmerer and Furbank, 1999), we assume, the rate of decarboxylation is equal to the rate of PEP carboxylation, so the same symbol, V * p , is used.
The diffusion of bicarbonate ions in epidermis, mesophyll cells, vasculature, or bundle sheath cells is collectively given by (Retta et al., 2016a): where D HCO3 − is the diffusivity of bicarbonate in the respective cell compartments. B CA and V * p is applied only to the mesophyll cytosol (Eqs 2, 3).
The steady state R-D equation for O 2 transport is given by: is the volumetric rate of oxygen production by chloroplasts of mesophyll or bundle sheath cells and zero elsewhere; and r * is the total O 2 consumption due to RuBP oxygenation during photorespiration in the bundle sheath chloroplast, and additional consumption of O 2 (0.5 mol of O 2 per RuBP oxygenation) by glycolate oxidase in peroxisomes in the photorespiratory cycle (Eq. 12) (Somerville, 2001). The rate of oxygen evolution, E * O2 , associated with the linear electron transport in chloroplasts was set to be equal to a quarter of the rate of linear electron flux (Eq. 19), because the production of 1 mol of O 2 is accompanied by the flow of four electrons.

Calculation of microscale model variables
Respiratory CO 2 release is assumed to occur in the epidermis, cytosol of mesophyll and bundle sheath cells, and the phloem tissue in the vascular bundles. Rate of respiratory CO 2 release in each of these components is calculated by assuming that the epidermis is as metabolically active as mesophyll, bundle sheath, or vascular bundles. Volumetric rate of respiratory CO 2 release in the epidermis (EPI), R EPI * , is given by: where R d is the rate of day respiration; f resp is the sum of volume fractions of vascular bundles, EPI, and cytosol of mesophyll and bundle sheath cells; and d is the average thickness of leaf tissue. Volumetric rate of respiratory CO 2 release in component i (mesophyll, bundle sheath, or vascular bundles), R i * , is given by: where f i is the volume fraction of a component i in the leaf tissue where the rate is applicable. The volumetric rate of photorespiratory CO 2 release in bundle sheath cytosol, r * p,CO2 , is given by: where V * c is given by Eq. 1; V BS,ch is the volume of bundle sheath chloroplast; and V BS,cy is the volume of the cytosol of bundle sheath. Volume of bundle sheath cytosol is calculated by ´V BS,cy (·)dv =˝dxdydz The average rate of CO 2 production from decarboxylation of C 4 acids, assumed to proceed at the same rate as PEP carboxylation (Eq. 2) in mesophyll cytosol (subscript M,cy) (von Caemmerer and Furbank, 1999), was calculated as: where V BS,ch is the volume of bundle sheath chloroplast (subscript BS,ch) and V M,cy is the volume of mesophyll cytosol. Consumption of 1 mol of O 2 by photorespiration results in the release of 0.5 mol of CO 2 . In addition, 0.5 mol O 2 is consumed per RuBP oxygenation by glycolate oxidase in peroxisomes in the photorespiratory cycle. Thus, the volumetric rate of oxygen consumption by photorespiration, r * p,O2 , in a bundle sheath chloroplast is given by: where the volume of bundle sheath chloroplast is calculated by ´V BS,ch (·)dv =˝dxdydz. Volumetric rates (marked with asterisk) of leaf area-based quantity Q such as V c,max or V p,max was calculated as: where f i is the volume fraction of a component i (mesophyll cytosol or bundle sheath chloroplast) in the leaf tissue where the rate is applicable. The rate of net photosynthesis was calculated by integrating the flux of CO 2 over the outer epidermis surfaces. The computed photosynthesis, A n, is given as: where Flux is a normal diffusive flux defined as n(−D∆C); n is the normal vector of unit length; S EPI is the epidermis surface; and S leaf is a projected area of the leaf. The rate of CO 2 leakage, L is given by: where V p is the mean rate of PEP carboxylation calculated from Eq. 2 and R M is the rate of CO 2 release in mesophyll cytosol through respiration (Eq. 9). The leakiness Φ is calculated as:

Calculation of fluxes and resistances to transport in boundaries
The flux of CO 2 , HCO 3 − , and O 2 at the boundary of cuticle, mesophyll, or bundle sheath cells is given by: where y stands for CO 2 , O 2 , and HCO 3 − ; P y is the permeability of cuticle, plasma membrane, chloroplast envelope, cell wall (CO 2 and O 2 ) or the mesophyll-bundle sheath interface to y; and ∆ [y] is the difference in concentration of y across each boundary. P y is given by: where D y is the diffusion coefficient of y in a cell component; ζ is the effective porosity; and t is the thickness of diffusion path. For the plasmodesmata, D y is assumed equal to that in the cytosol; ζ is the ratio of the bundle sheath surface area covered by plasmodesmata to the total surface area of bundle sheath; and t is the length of the plasmodesmata.
At the outer tangential cell wall of bundle sheath cells of maize, a suberin lamella is deposited (Hattersley and Browning, 1981). We assume that gas exchange occurs only through plasmodesmata at the mesophyll-bundle sheath interface. Thus, an insulated boundary condition was assumed at the surface of the bundle sheath cells exposed to the IAS.

ATP and NADPH production model
The ATP production was modeled to be driven by local light intensity at the level of the chloroplast using the light profile of the 3D geometry. The model of ATP production is explained in detail by Yin and Struik (2018). The model for electron transport and ATP production aimed to ensure that the ATP and NADPH requirements of C 4 and C 3 cycles were met at the whole-leaf level. A brief description of the model with selected equations relevant for this study is given below. The volumetric rate of linear electron transport (LET) in mesophyll or bundle sheath chloroplasts, j * i is given by a commonly used nonrectangular hyperbolic model (Farquhar and Wong, 1984;Yin et al., 2006): where i refers to either mesophyll or bundle sheath chloroplasts; α 2,LL,i is the photochemical efficiency of linear electron transport of photosystem II (PSII) under limiting light (LL) (based on light absorbed by both photosystems) when cyclic electron transport (CET) occurs simultaneously; I * abs,i is the rate of photon absorption by mesophyll or bundle sheath chloroplasts per unit leaf volume, calculated by the product of actinic irradiance (I inc ) and the fraction of that irradiance absorbed by the chloroplasts in these cells expressed per volume; j * max,i is the maximum volumetric electron transport rate under saturating illumination (calculated using Eq. 13); and θ is the convexity index. α 2,LL,i is calculated as (Yin et al., 2006): where φ 1,LL and φ 2,LL are the electron transport efficiencies under limiting light for photosystem I (PSI) and PSII, respectively (unitless, on the basis of light absorbed by each photosystem). f CET,i is the fraction of CET in mesophyll or bundle sheath chloroplasts given by: where j * LL,CET, i is the light-limited rate of CET, and j * LL,LET, i is the lightlimited rate of LET in either mesophyll or bundle sheath. j * LL,CET, i and j * LL,LET, i are calculated as : where u i is the fraction of light for LET in mesophyll or bundle sheath chloroplasts, which were calculated according to (Yin and Struik, 2018). The total ATP production is the sum of ATP production in mesophyll and in bundle sheath chloroplasts. The ATP production rate per leaf area basis in mesophyll or bundle sheath chloroplasts, j ATP,i, can be formulated as a function of PSII LET flux as : where h is the number of H + ions required to synthesize 1 ATP; f Q is the fraction of the electron flux that follows the Q cycle ; and S leaf is the area of a cross section of the computational domain. So, the term 1 Sleafí j * i dv represents the volumetrically integrated rate of LET per leaf area basis.
The rate of whole-leaf total ATP production calculated from chlorophyll fluorescence (CF) data, j ATP,CF , was given by (Yin et al., 2011): where sʹ is the lumped calibration factor; F s is the steady-state relative fluorescence yield; Fʹ m is the maximum relative fluorescence yield in the leaf; and x is the fraction of ATP allocated to the C 4 cycle. For this calculation, a constant x of 0.40, which is shown to be an optimal distribution for high photosynthesis, for all irradiances was assumed (von Caemmerer and Furbank, 1999;Yin et al., 2011;Yin and Struik, 2012). This value changes at low light intensity (Kromdijk et al., 2010). However, accounting for such change requires a detailed description of the kinetics of photosynthesis and exchange of metabolites, which is beyond the scope of our study (Bellasio, 2017).
The potential ATP production rate ( j ATP ) computed for the whole leaf level is given by: where S leaf is the area of a cross section of the computational domain; j * ATP is the local volumetric rate of ATP production determined from Eq. 24; and V leaf is the total leaf volume.
The light saturated electron transport rates in mesophyll and bundle sheath cells, j max,M and j max,BS (Eq. 19), were fitted by minimizing the difference between the respective ATP production rate determined from chlorophyll fluorescence measurement (Eq. 25) and that modeled from the potential electron transport rate at the whole leaf level (Eq. 26). To avoid overfitting, it was assumed that: where α 2LL,BS or α 2LL,M is given by Eq. 20; I abs,BS is the absorbed irradiance by bundle sheath chloroplasts; and I abs,M is the absorbed irradiance by mesophyll chloroplasts. The rate of NADPH production per leaf area in mesophyll or bundle sheath chloroplasts, j NADPH,i , is given by: where 0.5 refers to 0.5 mol NADPH produced per mol electrons transferred by LET. The fraction ATP in bundle sheath cells ( f ATP,BS) is calculated as: where j ATP,BS and j ATP,M are given by Eq. 24. The fraction of NADPH in bundle sheath cells ( f NADPH,BS) is calculated as: where j NADPH,BS and j NADPH,M are given by Eq. 28.

Modeling chloroplast arrangement
The chloroplast arrangements were modeled based on microscopic images of chloroplast movement for maize reported previously (Yamada et al., 2009;Maai et al., 2020b) assuming the vacuole allows free movement of chloroplasts. From the default geometry ( Supplementary Fig. S1C), two leaf geometries were developed to simulate the aggregative movement and the avoidance movement of mesophyll chloroplasts ( Supplementary Fig. S2D, E). The key morphological operations in MATLAB software (R2017b) used to model chloroplast arrangements were 'dilation', which adds pixels to the boundary of a domain, and 'erosion', which removes pixels from the boundary of a domain. The domain of bundle sheath cells of the default anatomy ( Supplementary Fig. S2C) was dilated to intersect the domain of mesophyll cells in contact with the bundle sheath cells to model the aggregative movement of mesophyll chloroplasts ( Supplementary Fig. S2D). The domain created as the result of the intersection was considered as the chloroplast domain ( Supplementary Fig. S2D) after it was eroded in MATLAB and the volume fraction of chloroplasts per cell matched reasonably well with that of the default geometry ( Supplementary Fig. S3). To model the avoidance movement of mesophyll chloroplasts ( Supplementary Fig.  S2E), a chloroplast layer was similarly made but only for those mesophyll cells lying between veins.

3D leaf anatomy, chloroplast arrangement, and ATP and NADPH production
The 3D model of light propagation was coupled to the model of ATP and NADPH production (Eqs 19-24, 28) to determine the total ATP and NADPH production. The total ATP was partitioned between C 3 and C 4 cells at cell level. Forty percent of the total produced ATP (Eq. 24) was assumed to be used for the C 4 cycle (von Caemmerer and Furbank, 1999). Then, the gas transport models (Eqs 5-7) were solved using the 3D structure for various chloroplast arrangement scenarios. The responses of A n (Eq. 14) to I inc and to intercellular CO 2 (C i ) were computed and compared with those obtained from the gas exchange measurement. The fractions of CET in bundle sheath cells (Eq. 21) and the fraction of ATP and NADPH production in bundle sheath chloroplasts (Eqs 29-30) calculated from the default geometry were compared with cases of chloroplast arrangement.

Comparison of a 2D and a 3D model of gas transport
Three 2D slices from the 3D geometry were taken from the positions 30, 60, and 90 of 120 slices to directly compare the differences in modeled photosynthesis on the same leaf using the same R-D approach between the simpler 2D geometry and the corresponding 3D geometry ( Supplementary Fig. S4). The light absorption profile and the mean ATP production rate were set equivalent to those of the 3D geometry at the positions where the slices were made. Following Retta et al. (2016a), we assumed that the concentration of CO 2 in the IAS was uniform for the 2D model. Therefore, the gas concentrations at the exposed surface of mesophyll were set equal to the mean concentrations in IAS of the 3D model. Retta et al. (2016a) also assumed that the diffusion coefficients of CO 2 and bicarbonate were equal to those of pure water in order to improve the match between the measured and simulated gas exchange data. Thus, we also solved the 2D model from slice 60 with these assumptions so that the relative viscosity of the cell cytosol (η) was assumed to be 1.0 or 2.0, where the latter value is assumed closest to the biological medium. All other parameter values were kept the same as the ones reported in Table 1. The response of A n to C i was computed for the 2D model (for various assumption of η) and compared with the results obtained from the 3D model.

Numerical solution
A 3D leaf tissue of 124 × 124 × 200 µm was represented by 9.428× 10 6 cube elements with a length of 0.67 µm. The lateral sides of this geometry were assumed to have no net flux, while external concentrations of CO 2 and O 2 were assumed at the adaxial and abaxial sides. The R-D equations were discretized over the finite volume grid and the resulting linear systems of algebraic equations of the unknown concentrations at the volume nodes were solved by the preconditioned conjugate gradient procedure in MATLAB using the finite volume method (Versteeg and Malalasekera, 1995;Ho et al., 2016). We used a 16-GB RAM node of the high-performance computer in the VSC-Flemish Supercomputer Center (Belgium) to run the model. The computation time was 72 h.

Results
The 3D leaf structure impacts gas transport and light absorbance The 3D geometry (Fig. 1A) consisted of the epidermis, stomata, IAS, mesophyll cells, bundle sheath cells and vasculature, as main components of the maize leaf anatomy. Leaf anatomical characteristics measured from this geometry were compared with literature values ( Table 2). The anatomical features such as exposed surface area of mesophyll per leaf area, surface area of bundle sheath cells per leaf area, surface to volume ratio of mesophyll and bundle sheath cells and volume fraction of chloroplasts in mesophyll and bundle sheath compared well with literature. Table 2 shows that the volume per leaf area of mesophyll and bundle sheath cells of the 3D geometry were lower than those reported in the literature. Figure 2A shows the comparison of optical properties measured experimentally with those predicted from the light propagation model using the default mesophyll chloroplasts configuration. The root mean square errors for the prediction of spectra of optical properties were 4% and 10% for total reflectance and transmittance ( Fig. 2A), respectively. The simulation model for light propagation thus gave reasonable predictions for the reflectance and transmittance of red light (which was 90% of the total applied light in the gas exchange measurement, Fig. 2A). The fraction of photon energy absorbed by the leaf (Fig. 2B, C), as illustrated by a 3D representation in Fig.  1D, shows that halfway through the leaf depth, the photon absorption was enhanced by the chloroplasts within the bundle sheath cells. There was little effect of changing the scattering coefficient of chloroplasts by 25% on the absorption profile suggesting the effectiveness of our approach of accounting the scattering property of cell wall with that of chloroplasts. The response of j ATP to I inc (Fig. 2D) was computed using the estimated maximum rate of electron transport (Table 1). The light model predicts that 87% was absorbed by the leaf, while the mesophyll chloroplasts absorbed 55% and the bundle sheath chloroplasts absorbed 28% (Table 3). The predicted responses of A n to changes in I inc and C i are shown in Fig. 3A, B. The mean square percentage error for the A n -I inc curve was 6.5% and that of the A n -C i curve was 1.8%. The deviations at low light intensity contributed to the higher mean square percentage error of the A n -I inc curve. At high light intensity where the chloroplast movement scenarios were considered, the error in the prediction was much lower (Fig. 3). Thus, the model could reasonably predict the response curves. The mean [CO 2 ] in the bundle sheath chloroplast (C c ) was about five times higher than ambient [CO 2 ] (380 µmol mol −1 ) when computed at 210 mmol mol −1 O 2 and I inc of 1500 µmol m −2 s −1 due to the CCM (Fig. 3C). The predicted CO 2 concentration profile in the IAS was largely homogeneous (Fig. 3D).
The 3D model was compared with a simpler 2D model (Retta et al., 2016a) using random slices from the same default leaf geometry (Supplementary Fig. S4). The A n -C i curve was underestimated by the 2D models, especially if the viscosity of the aqueous symplast was estimated as double that of pure water (Supplementary Fig. S4A). The C c estimated from the 3D model was 1584 µmol mol −1 (ambient [CO 2 ] of 380 µmol mol −1 , 210 mmol mol −1 O 2 , and I inc of 1500 µmol m −2 s −1 ). The C c calculated using three 2D slices selected at z-positions of 20, 40, and 60 in the stack of 120 slices was 1520 ± 270 µmol mol −1 suggesting that if C c estimated from the 3D model is the more accurate, the use of randomly selected 2D slices will bias the concentration of CO 2 at the sites of Rubisco. In particular, the 3D model of light absorption within the leaf structure profile (Fig. 1D) represents a major improvement over the assumption of constant photon flux across the leaf in the model of Retta et al. (2016a). Moreover, 2D models understate the interconnectivity of the IAS prevalent in the 3D space and may lead to areas of IAS not connected to a stoma, which compromises the modeling of CO 2 diffusion within the IAS (Supplementary Fig. S4). Therefore, the 3D model is definitely more accurate than the 2D model to simulate diffusion of gases and photosynthesis in a C 4 leaf.

Chloroplast movement increased the relative light absorbance by bundle sheath chloroplasts
wThe effect of mesophyll and bundle sheath chloroplast arrangement on light propagation is illustrated in Fig. 1D-F. The geometries of the different chloroplast arrangement scenarios (Fig. 1B, C) have similar volume fractions of chloroplasts per cell and total volume fractions of chloroplasts as the default geometry ( Fig. 1A; Supplementary Fig. S3). The interveinal mesophyll cells whose chloroplast arrangement changed accounted for 25% of the mesophyll volume and the total chloroplast in avoidance arrangement is 18% of the total chloroplast volume. Table 3 shows the total light absorption by the leaf decreased by ca. 7% for aggregative movement of mesophyll chloroplasts and ca. 6% for avoidance movement compared with that of the default geometry. The avoidance movements reduced the absorbance of chloroplasts by ca. 15% while aggregative movement resulted in ca. 8% lower absorbance. On the other hand, the transmittance of light through the maize leaf doubled in both chloroplast movement scenarios (Table 3). By contrast, the ratio of absorbed light by bundle sheath chloroplasts relative to mesophyll chloroplasts was increased by ca. 73% over the default case due to the aggregative movement and ca. 45% due to avoidance movement. Although transmittance did not change (Table 3) the reflectance increased by 6% in the avoidance movement compared with the aggregative movement.
The total leaf absorbance profile across the leaf depth is shown in Fig. 2B. The profile for the avoidance movement slightly differed close to mid-leaf thickness from that of the default case. The profile for avoidance movement changed in such a way that there was an increased absorbance followed by steeper decline starting at 100 µm depth, which was where the chloroplast arrangement changed from that of the default geometry (Fig. 2B). The aggregative chloroplast arrangement had a more significant change in the profile of light absorbed across the leaf, which resulted in a peak of absorption around the bundle sheath (Fig. 2B), but also boosted the fractional photon absorption by bundle sheath chloroplasts (Fig. 2C). Although the fractional photon absorption by bundle sheath increased in the avoidance scenario at lower leaf deep, it decreased deeper in the leaf (Fig. 2C).

Aggregative movement increased ATP and NADPH production in bundle sheath at high light intensities
The effect of chloroplast arrangement on ATP and NADPH production required to balance cell-type-specific demand at limiting light conditions is shown in Table 3. For the default geometry, the absorbed light was mainly used to drive linear electron transport (LET) in mesophyll chloroplasts and cyclic electron transport (CET) in bundle sheath cells. About 40% of the total ATP production was due to bundle sheath chloroplasts. The low NADPH production in bundle sheath cells was due to the assumption of little PSII in bundle sheath cells in the NADP-ME species.
In the aggregative movement case, nearly all of the light absorbed by mesophyll chloroplasts should be used for LET, as in the default case (Table 3). However, ca. 18% more light was allocated to LET in bundle sheath chloroplasts compared with the default. At the bundle sheath cell level, the contribution of CET to the total electron transport rate was only slightly lower in the aggregative than in the default arrangement (Table  3). The aggregative arrangement required ca. 36% more ATP production and ca. 64% more NADPH production in bundle sheath cells.
For avoidance movement, there was increased use of absorbed light by mesophyll and bundle sheath chloroplasts for LET. The increased use of light for LET was ca. 8% more in mesophyll chloroplasts and ca. 11% more in bundle sheath chloroplasts. In addition, the fractions of ATP and NADPH production at limited light intensity (Table 3) increased by ca. 43% and ca. 26% respectively, compared with the default chloroplast arrangement.
Aggregative chloroplast movement increased A n at high light but also increased leakiness Figure 4A shows the effect of chloroplast arrangement on light response of A n . Because of the higher ATP production (Fig. 2D), aggregative movement resulted in 4-8% higher A n at high I inc where the movement may occur. By contrast, A n decreased by ca. 12-16% for avoidance movement of mesophyll chloroplasts at I inc of 500 µmol m −2 s −1 and higher. At I inc of 1500 µmol m −2 s −1 , Φ was ca. 30% higher for the aggregative and avoidance arrangements than that for the default geometry (Fig. 4B). The increase in Φ for aggregative movement was due to higher CO 2 in bundle sheath cells and leakage created by the imbalance in CO 2 assimilation in mesophyll and bundle sheath ( Supplementary Fig. S5). The higher ATP production at high I inc translated into increased PEP carboxylation as the latter was limited by PEP regeneration, which consumes ATP. The higher PEP carboxylation increased the CO 2 supply into bundle sheath chloroplasts, where the local Rubisco carboxylation was not proportionally increased ( Supplementary Fig. S5). In addition, aggregative and avoidance movements increased oxygen production in bundle sheath cells because of the higher LET (Table 1) thereby increasing photorespiratory CO 2 release in bundle sheath cytosol (not shown) despite a lower PSII proportion in bundle sheath chloroplasts of maize (assumed here to be 10%).

Overall performances of the 3D reaction-diffusion model for C 4 photosynthesis applied to a maize leaf
We developed and tested for the first time a model that concurrently simulates CO 2 and O 2 transport, light propagation, C 4 energetics and related photosynthetic processes in the full 3D geometry of a leaf of the C 4 model crop maize. The 3D model was primarily validated by comparing simulated and measured responses of photosynthesis to CO 2 concentration and light intensity. A large part of the merit in this rests on the fine tuning of the 3D light propagation model within a real maize leaf that properly reproduced the in vivo absorptance and transmittance measurements. This first 3D model of C 4 photosynthesis constitutes an advance on the understanding of how anatomy regulates photosynthesis in C 4 leaves, a complex process that requires fine tuning between the CCM and light available for chloroplasts in mesophyll and bundle sheath cells. The anatomical features influencing gas diffusion and light propagation compared well with literature values, although the features vary with treatments (Table 1) (Dengler et al., 1994;von Caemmerer et al., 2007;Retta et al., 2016b). The 3D model was applied in analysing photosynthesis by addressing the role of chloroplast arrangement, which requires a representative tissue-level 3D geometry and could not be realistically carried out by a 2D model ( Supplementary Fig. S4) (Théroux-Rancourt et al., 2017;Earles et al., 2019). The 3D light profile showed that large photon absorption gradients existed when the Kranz anatomy was illuminated unilaterally depending on the 3D chloroplast distribution across the leaf thickness.
The large gradient means that some of the chloroplasts are still light limited (Figs 1D, 2B; 18% by volume of the chloroplast absorbed less than 0.001% of the photons) and may strongly contribute to the non-saturated irradiance response of the leaflevel photosynthesis. However, light may be scattered more than modeled here by cell walls (Gates et al., 1965;Allen et al., 1973;Vogelmann, 1993), which were not explicitly included in the geometry because of the requirement of a much denser mesh and, hence, computational demand, but the scattering effect was lumped with the chloroplasts. Light may also penetrate deeper in leaf tissue in C 4 leaves having bundle sheath extensions (Bellasio and Lundgren, 2016), which were not included here because bundle sheath extensions were not observed in the maize leaf anatomy (Supplementary Fig. S1A). In addition, the 3D model could be validated more accurately by comparing the intra-leaf profile of light and photosynthetic capacity. For instance, light-sheet microscopy can be used to resolve the profile of light across the leaf thickness (Slattery et al., 2016). In addition, the profile of photosynthetic capacity could also be estimated using epi-illumination fluorescence microscopy to derive the profile of chlorophyll and Rubisco, together with the aforementioned light absorption profile (Vogelmann and Evans, 2002;Slattery et al., 2016;Borsuk and Brodersen, 2019). Overall, the heterogeneity in light and [CO 2 ] suggests that the contributions of electron transport and enzyme-limited CO 2 assimilation rates to the total assimilation are influenced by the leaf microstructure.
Effect of mesophyll chloroplast movement on CO 2 concentrating mechanism efficiency and photosynthetic rate The 3D gas transport model was coupled with the ATP and NADPH production model using the validated light propagation model and was applied to test hypotheses on the effect of chloroplast movement on CCM efficiency and photosynthetic rate. The default 3D geometry was modified to model the geometries of the various arrangements of mesophyll and bundle sheath chloroplast scenarios inspired by what is observed in vivo from microscopy imaging (Yamada et al., 2009;Maai et al., 2011Maai et al., , 2020aSales et al., 2018). Care was taken not to modify the chloroplast volume fraction at the cell level in addition to the total chloroplast volume fraction so that the absorption and scattering coefficients remained similar. Consistent with previous observations, chloroplast movement increased transmittance and decreased light absorption  (Yamada et al., 2009;Sun et al., 2014;Maai et al., 2020a). Avoidance movement was slightly more effective in decreasing total light absorption by chloroplasts. By contrast, aggregative movement highly increased absorption by bundle sheath chloroplasts while reducing the light exposure of mesophyll chloroplasts (Table 3). Although the aim of this study was to test the efficiency of mesophyll chloroplasts in gas exchange rates, our findings of light absorption profiles could be used in further studies to evaluate mechanisms of photoprotection and photoinhibition in C 4 leaves when abiotic stress results in chloroplast movement (Ghannoum, 2009;Guidi et al., 2019). The 3D model predicts that chloroplast movement impacted photosynthetic rate and CCM efficiency through altered light availability in mesophyll and bundle sheath chloroplasts, energy production, and changes in the process of electron transport. At high light intensities, where the aggregative movement was expected, there was a strong increase in photosynthetic rate (Fig. 4A) supporting our hypothesis that aggregative movement of mesophyll chloroplasts improves photosynthetic rate. In addition, avoidance movement of mesophyll chloroplasts decreased ATP production and A n and increased Φ at high light intensities compared with the default case of no chloroplast movement. Thus, the hypothesis that avoidance movement decreases photosynthetic rate was supported. These results support the previous suggestion that preference for aggregative movement rather than avoidance movement in maize upon exposure to blue light intensity may have benefits for photosynthesis (Ryu et al., 2014). The altered light absorbance by mesophyll and bundle sheath chloroplasts created an imbalance in the rate of CO 2 assimilation in mesophyll and bundle sheath cells resulting in increased leakiness (Fig. 4b). Such an imbalance has recently been shown to increase leakiness due to lack of coordination between mesophyll and bundle sheath cells following dark to high light transition (Wang et al., 2022). Previous modeling results suggested that increased leakiness may be beneficial as it helps balance reducing power demand and availability by removing excess CO 2 from bundle sheath cells (Wang et al., 2014;Kromdijk et al., 2014). Balancing ATP and NADPH demand of the CCM required an increased operation of linear electron transport in mesophyll cells due to less light absorption by mesophyll chloroplasts (Table 3). Overall, the results support the hypothesis that chloroplast movement, although useful for photoprotection because less light is absorbed by chloroplasts, may affect the operation of CCM and reduce its efficiency (Taniguchi et al., 2003;Yamada et al., 2009;Maai et al., 2011Maai et al., , 2020bSun et al., 2014). These results are clear indications that C 4 photosynthesis can be improved by optimizing light absorption and increasing ATP and NADPH production (Ermakova et al., 2019(Ermakova et al., , 2021b, but at the expense of CCM efficiency. These results are relevant for stacking traits in the C 4 Rice project (Ermakova et al., 2021a;Xiao et al., 2023), for example, one should consider designing strategies of boosting Rubisco activity in order to remedy the lower CCM efficiency in addition to reducing the limitation by electron transport at high light intensity by boosting electron transport capacity in C 4 plants (Ermakova et al., 2019).
Chloroplast movement in C 4 plants, unlike in C 3 plants, requires exposure to high light intensity for several hours, which can be experienced by leaves close to the top of the canopy. The long duration may have been due to the sparse arrangement of chloroplasts in mesophyll cells of C 4 plants making them less susceptible (Ghannoum et al., 2005;Stata et al., 2014), an evolutionary adaption of C 4 plants, which may allow more light to pass through to bundle sheath cells (Stata et al., 2014). Our results suggest that CCM efficiency may also add another constraint, which needs to be verified experimentally. In addition, the limitation of metabolite transport may explain the limited movement of maize bundle sheath chloroplasts (Taniguchi et al., 2003;Kato et al., 2022). Our simulation results suggest that aggregative movement may provide benefits for improved photosynthesis besides photoprotection. However, if mesophyll chloroplasts are not arranged in this mode most of the day, it means that random distribution of chloroplasts in mesophyll cells is needed for processes other than photosynthesis, and only when light intensity is really high does the advantage of this configuration for gas exchange and photoprotection overcome the benefits of the other chloroplast functions within mesophyll cells.

Concluding remarks
We developed a 3D R-D model of gas exchange and light propagation coupled to an ATP and NADPH production model and applied it to test the hypotheses that movement of chloroplasts reduces light absorption within the C 4 leaf but modifies differently the CCM efficiency and photosynthesis depending on the type of mesophyll chloroplast movement. We found that mesophyll aggregative chloroplasts arrangement increased the proportion of light absorption by bundle sheath cells, the leaf ATP and NADPH production, and photosynthetic rate. By contrast, the avoidance arrangement decreased photosynthesis. None of the chloroplast movement scenarios were beneficial for maintaining CCM efficiency. C 4 plants could boost photosynthetic rate further at the expense of increased leakiness using chloroplast movement in addition to photo-protection under high light intensity (Yamada et al., 2009;Maai et al., 2011Maai et al., , 2020b. Furthermore, evidence has been provided that a 3D model, unlike simpler biochemical models and a 2D model, can provide increased mechanistic understanding of the light propagation, CO 2 diffusion, and chloroplast movement in complex C 4 photosynthesis.

Supplementary data
The following supplementary data are available at JXB online. Fig. S1. Sample light microscopy images of maize leaf tissue. Fig. S2. Stomatal conductance to CO 2 in response to irradiance and degree of stomatal opening.   Fig. S4. Comparison of a 2D and a 3D model. Fig. S5. Comparison of local rates and CO 2 concentration. Table S1. List of symbols, their definitions, and units. Table S2. Mean equivalent radius and total number of organelles per volume. Table S3. Computed optical properties of the different compartments of the leaf model.