In silico evidence for the utility of parsimonious root phenotypes for improved vegetative growth and carbon sequestration under drought

Drought is a primary constraint to crop yields and climate change is expected to increase the frequency and severity of drought stress in the future. It has been hypothesized that crops can be made more resistant to drought and better able to sequester atmospheric carbon in the soil by selecting appropriate root phenotypes. We introduce OpenSimRoot_v2, an upgraded version of the functional-structural plant/soil model OpenSimRoot, and use it to test the utility of a maize root phenotype with fewer and steeper axial roots, reduced lateral root branching density, and more aerenchyma formation (i.e. the ‘Steep, Cheap, and Deep’ (SCD) ideotype) and different combinations of underlying SCD root phene states under rainfed and drought conditions in three distinct maize growing pedoclimatic environments in the USA, Nigeria, and Mexico. In all environments where plants are subjected to drought stress the SCD ideotype as well as several intermediate phenotypes lead to greater shoot biomass after 42 days. As an additional advantage, the amount of carbon deposited below 50 cm in the soil is twice as great for the SCD phenotype as for the reference phenotype in 5 out of 6 simulated environments. We conclude that crop growth and deep soil carbon deposition can be improved by breeding maize plants with fewer axial roots, reduced lateral root branching density, and more aerenchyma formation.


Introduction
Globally, the incidence and severity of droughts under climate change have been increasing over the past few decades, causing substantial agricultural losses with cascading enviro-socioeconomic impacts (Reddy and Hodges, 2000;Mukherjee et al., 2018). Drought is a complex phenomenon commonly characterized by suboptimal availability of water for a sustained period (Wilhite and Glantz, 1985;Van Loon, 2015). In an agricultural context, drought is attributed to the reduction in water supply required by the plant due to soil water deficit, mostly in the root zone, and/or an increase in water loss via transpiration (Boken et al., 2005;Salehi-Lisar and Bakhshayeshan-Agdam, 2016). Soil water availability and plant adaptation to water deficit are spatiotemporally influenced by an array of pedoclimatic factors. These factors pose a major challenge in understanding drought and its impact on plants, which has implications for improving agricultural practices and breeding efforts.
Plant water status reflects the synchronized response of the plant to soil water availability and atmospheric demand (Silva and Lambers, 2021). Shoot water loss is driven by atmospheric demand together with leaf area, stomatal conductance, and intrinsic water uptake capacity of the plant (De Swaef et al., 2022). On the other hand, soil properties and root architecture largely determine the water availability to the plant. To avoid water deficit, plants elicit an array of response mechanisms, which broadly involve -a) limiting water loss via stomatal adjustment and limiting shoot growth, b) increasing root foraging to enhance water uptake, and c) restricting tissue dehydration via osmotic and metabolic adjustments (Bray, 1997;Dodd and Ryan, 2016;Salehi-Lisar and Bakhshayeshan-Agdam, 2016).
One avenue to counter water loss via transpiration is to increase root water capture. Since water is often available in deeper soil domains, in most environments deeper rooting improves drought resistance. In this context, the 'Steep, Cheap, and Deep' (SCD) root ideotype was proposed for improving drought resistance in crops, by increasing root depth and in turn water acquisition from deep soil domains (Lynch, 2013). Fundamentally, the SCD ideotype is an integrated root phenotype that aims to optimize how internal plant resources and soil foraging activities are spatiotemporally distributed among different root classes and across the plant to improve its performance in response to nutrient stresses and drought. With the root system adapted at morphological, anatomical, and physiological scales, the SCD ideotypes that enable root foraging of deeper soil strata are typically characterized by steeper axial root angles, reduced number of axial roots, reduced lateral branching density, and formation of root cortical aerenchyma and other anatomical phenotypes that reduce the metabolic cost of soil exploration.
Over the years, various other combinations of root phenotypes corresponding to SCD ideotypes have been proposed that, in silico and/or in vivo, have been associated with improved crop growth under suboptimal water and nitrogen availability, in various species, such as barley, maize and rice (Lynch 2018;Lynch, 2019;Lynch, 2022). Reducing the number of axial roots decreases competition among axial roots for internal plant resources and external soil resources, which allows the remaining axial roots to grow deeper, increasing water capture from deeper soil layers. This was confirmed in a field study where maize plants with fewer axial roots had increased rooting depth, leaf water content, shoot biomass, and yield under drought (Gao and Lynch, 2016). Nitrogen in agricultural soils is generally a mobile resource (as nitrate) that, like water, is often more available in deep soil strata (Thorup-Kristensen et al., 2020). Reduced axial root number increased shoot biomass for maize and rice in response to suboptimal nitrogen availability (Saengwilai et al., 2014;Gao and Lynch, 2016;Ajmera et al., 2022). When the soil gets very dry and soil hydraulic conductivity becomes extremely low, there is significant competition for water and nitrogen among neighboring lateral roots. This means that reducing lateral root branching density is beneficial for plant performance under suboptimal water and nitrogen availability . Furthermore, the utility of steeper axial root angles for improved tolerance to drought and low nitrogen supply by developing a deeper root system has been confirmed in different crop species (Manschadi et al., 2008;Trachsel et al., 2013;Uga et al., 2013;Lynch and Wojciechowski, 2015;Schneider et al., 2022). Since one of the secondary effects of water deprivation on plants is a reduction in photosynthesis, plants suffering from water shortages will have limited photosynthetically-derived carbon. Anatomical phenotypes that reduce the metabolic cost of soil exploration should therefore improve soil exploration and water capture under drought, which is the case with root cortical aerenchyma, reduced cortical cell file number, and increased cortical cell size in maize, rice and wheat (Jeong et al., 2013;Slack 2018;Lynch et al., 2021).
It is noteworthy that a substantial amount of photosynthate is allocated to roots and eventually deposited into the rhizosphere. Shoot-derived carbon decays more rapidly than root-derived carbon (Puget and Drinkwater, 2001;Ghafoor et al., 2017). In addition, oxygen availability and microbal activity both tend to decrease with soil depth, slowing the decomposition of plant residues. Given this, root phene states that enable greater rooting depth and thus increase soil carbon deposition could be useful in removing atmospheric CO 2 and concurrently, increasing soil organic matter (Kell, 2011(Kell, , 2012Lynch, 2015).
The SCD phenotype is composed of multiple root phene states. These phene states interact and influence the utility of each other. The interactions among different phene states vary with environments including water regimes and pedoclimatic conditions. Combining a substantial number of phenotypes with different environments creates a large number of different scenarios. Given this, simulation modeling is the only feasible way to explore the 'fitness landscape' of root phenotypes. Several functional-structural plant models have been developed which explicitly simulate root systems to explore and understand their interactions with soil properties and resources; and in turn, their impact on shoot growth. This includes models such as ArchiSimple (Pagès et al., 2014), CPlantBox (Zhou et al., 2020), OpenSimRoot , R-SWMS (Javaux et al., 2008) and SPACSYS (Wu et al., 2007). Besides root-soil interactions, the functional-structural plant/soil model OpenSimRoot simulates carbon assimilation, carbon partitioning for tissue growth and maintenance; and the effect of carbon restrictions on tissue growth and development . As a consequence of this capability to capture carbon economy, plant growth in OpenSimroot is determined both by the specified phenotype and its interactions with environmental conditions. In the past, OpenSimRoot has been useful in evaluating the utility of several root phenotypes for improving plant performance under low soil nitrogen, phosphorus, and potassium availability for different crop species in varying soil and climatic conditions (Postma and Lynch, 2011;Postma et al., 2014;Zhang et al., 2014;Dathe et al., 2016;Schneider et al., 2017;Rangarajan et al., 2018;Ajmera et al., 2022).
In this work, we introduce OpenSimRoot_v2, which adds new functionality to OpenSimRoot, enabling simulation of plant and soil responses to drought. This functionality includes the combined implementation of different bio-physio-chemical models available in the literature; namely the Farquhar-von Caemmerer-Berry (FvCB) model for C 3 and C 4 photosynthesis, leaf gas exchange and stomatal conductance (Von Caemmerer and Farquhar, 1981;Leuning, 1995), leaf temperature and energy balance models (Von Caemmerer, 2000), sun/shade model for leaf-to-canopy scaling (De Pury and Farquhar, 1997), a model for nitrogen-limited photosynthesis (Kull and Kruijt, 1998), water stress response functions, and models for simulating day-night cycles and the corresponding carbon allocations (Sulpice et al., 2014). With this upgrade, the OpenSimRoot_v2 still remains backward compatible enabling implementation of earlier OpenSimRoot versions with an appropriate selection of models and functions in the input files.
To better understand the relationship between root phenotypes and drought resistance, this in silico study evaluates the "Steep, Cheap, and Deep" ideotype and all combinations of the associated phene states in maize growing environments in the USA, Nigeria, and Mexico with distinct pedoclimatic conditions under rainfed and terminal drought conditions. Our objectives are: • Ascertain if the "Steep, Cheap, and Deep" ideotype is advantageous, that is, leads to greater shoot biomass, under water deficit stress (henceforth simply 'drought stress'). • Determine the extent to which interactions among phene states of the "Steep, Cheap, and Deep" ideotype contribute to performance under drought.
• Quantify the effects of different root phenotypes on shoot biomass, rooting depth, estimated carbon deposition at depth, and carbon use efficiency (water uptake per unit of carbon invested in roots) in rainfed and drought conditions.

Materials and methods
OpenSimRoot is a feature-rich, open-source, functionalstructural plant/soil simulation model, which explicitly simulates the geometry of roots and soil properties in three dimensions over time. In OpenSimRoot, the shoot is simulated through a number of state variables in a non-geometrical fashion. Being modular, OpenSimRoot encapsulates different submodels to represent various processes such as nutrient and water acquisition, aerenchyma formation, carbon sources and sinks, and nitrogen mineralization. Plant growth is determined by potential growth rates, constraints imposed by the availability of nutrients and carbon, and modifiers depending on the environment. These particular environments include soils differing in bulk density, texture, and hydraulic properties.
OpenSimRoot implements an ensemble of three models to capture water movement through the soil and the plant and into the atmosphere. This includes the SWMS3D model (Simunek et al., 1995) for soil water dynamics, a hydraulic network model (Alm et al., 1992;Doussan et al., 1998) for water uptake and transport through roots and the Penman-Monteith equations for evapotranspiration (Penman 1948;Monteith 1965). The SWMS3D model simulates water transport in the 3D soil by solving the Richards equation, which combines Darcy's law with mass conservation, using a finite element method . Soil water dynamics are determined by soil hydraulic properties, precipitation, evaporation, root water uptake as well as the specified boundary conditions at the bottom of the soil column. Soil bulk density and van Genuchten parameters determined by soil texture and organic content, define soil water the relationship between water content and matric potential for a given soil. These parameters vary with soil depth and are interpolated between the depths to avoid discontinuities along the border between two soil layers. Using the root impedance model added to OpenSimRoot in (Strock et al., 2022), we capture the impact of soil bulk density on root growth. Also included are models that consider the effect of soil water status on soil penetration resistance, with wetter soils being easier to penetrate. This is especially important when simulating plants in drying soil since root distribution affects the distribution of soil water and vice versa.
To accurately simulate the effects of drought stress on stomatal conductance, photosynthesis, and overall plant development we introduce here OpenSimRoot_v2, which adds submodels for various state variables related to photosynthesis. Figure 1 provides an overview of the most important submodels which were added and updated and their most important interactions. Note that many submodels and interactions were left out of this diagram for the sake of readability. Table 1 contains an alphabetically sorted list of symbols for derived variables related to the new or updated models with units and the relevant equation or section. Table 2 contains an alphabetically sorted list of symbols for derived quantities already previously present in OpenSimRoot with units. Table 3 contains an alphabetically sorted list of model inputs with the values used for this paper, units and references.
In OpenSimRoot, state variables are represented by objects which can be constant, follow a predefined variation over time, drawn from a random distribution, or calculated by a relevant class. Information moves between different parts of the model through a common application programming interface (API). Because of this modular structure, one can easily switch between different model implementations for a given state variable without the rest of the model being affected. Likewise, adding new state variables and model implementations is relatively straightforward and usually does not require any modification to the rest of the code. The new models added here required some additions to the OpenSimRoot engine, which we describe in section 2.9. We added new state variables and corresponding models and added updated models for some existing state variables to OpenSimRoot, which we describe in the following subsections. Readers are encouraged to consult the referenced papers for more details on and explanations of the model equations in these subsections.

Water stress
To properly model the effect of drought on plant functioning and development we need to quantify water stress. We do this by adding a water stress minimodel similar to how nutrient stresses are captured in OpenSimRoot. This water stress minimodel captures water stress as a number S w between 0 and 1, where 1 is no stress and 0 is maximum stress. It is known that abscisic acid (ABA) plays a role in the plant drought response (Kriedemann et al., 1972;Loveys, 1977;Ikegami et al., 2009;Okamoto et al., 2009;McAdam and Brodribb, 2016;Sack et al., 2018), but OpenSimRoot does not simulate the dynamics of phytohormones neither does it explicitly capture shoot geometry to simulate leaf water potential. As a result, we use a proxy to leaf water potential to capture drought stress in OpenSimRoot_v2. Notably, OpenSimRoot does includes a xylem water flow model described in (Doussan et al., 1998;Postma et al., 2017), which assumes steady-state flow to calculate the quantity and An overview of important new models and their dependencies in OpenSimRoot. Arrows point from a variable to variables depending on it, arrow colors are solely for distinguishing between them. White indicates models which were already present and are unchanged, red indicates models which were updated, blue indicates newly added models, green indicates newly added models with separate sunlit and shaded instances and purple indicates already present models which now have separate sunlit and shaded instances. distribution of water uptake by the roots for a given hydraulic potential at the collar (which is where the hypocotyl and shoot meet). This collar hydraulic potential is set such that the water uptake by the root system matches the potential transpiration, which is either calculated using the Penman-Monteith equations (Zotarelli et al., 2010) or based on photosynthesis rates. In OpenSimRoot_v2, we translate the collar potential Y c in hPa to a stress factor S w (Y c ) using a drought response curve defined in the input files. In this paper we used the following drought response curve: The precise relationship between leaf water potential and responses related to drought stress is difficult to establish but some estimates are available in the literature (Acevedo et al., 1971;Watts, 1974;Tanguilig et al., 1987). The water stress factor is translated to plant responses through transfer functions defined in the input file.

Photosynthesis
In order to have an accurate plant response to drought, the photosynthesis model we use should depend on solar irradiation as well as stomatal conductance. The Farquhar-Von Caemmerer-Berry model Caemmerer (Farquhar et al., 1980;Von Caemmerer and Farquhar, 1981;Farquhar and Caemmerer, 1982) satisfies these conditions, is the most widely cited photosynthesis model and many published steady-state photosynthesis models are based on it or derived from it (Morales et al., 2018). This model combines biochemical and irradiation dependent constraints and has been applied in a wide range of contexts (Von Caemmerer, 2013). For C 3 plants the photosynthetic assimilation rate A in mmol m 2 s is equal to where A c is the RuBisCo-limited (or carbon-limited) rate, A j the electron transport-limited (or light-limited) rate and A p the phosphate-limited rate of carbon assimilation, all in mmol m 2 s . We will assume that phosphate is not limiting so ignore A p from here onward. A c and A j are defined as Here V cmax is the maximum RuBisCo carboxylation rate in mmol m 2 s , K c and K O are the RuBisCo Michaelis constants for CO 2 and O 2 in mmol mol and Г* is the CO 2 compensation point without dark respiration (the leaf respiration that happens independently of photosynthesis) in mmol mol . C m and O m are the mesophyll CO 2 and O 2 concentrations in mmol mol (see Section 2.3). J is the potential electron transport rate in mmol m 2 s , which is defined as: where I 2 is the amount of energy in the light absorbed by photosystem II in mmol m 2 s , J max is the maximum electron transport rate in mmol m 2 s and q is an empirical curvature factor. I 2 is given by where b is a conversion factor in mmol j I s is the solar irradiation in W m 2 , a the leaf absorptance and f a factor to correct for the spectral quality of the light. The above equations are for C 3 photosynthesis, we use the following expressions for C 4 light-limited and carbon-limited photosynthesis rates (Von Caemmerer, 2000): Here C s and O s are the bundle sheath carbon dioxide and oxygen concentrations in mmol mol , S c/o is the RuBisCo specificity for CO 2 relative to O 2 and x is the electron transport rate partitioning factor.

Leaf gas concentrations
Photosynthesis rates depend on the CO 2 and O 2 concentrations in mesophyll (C 3 photosynthesis) or bundle sheath (C 4 photosynthesis) cells. We formulate simple models for these concentrations. In the case of C 3 photosynthesis, we model the mesophyll CO 2 and O 2 concentrations, C m and O m in mmol mol , using ordinary differential equations, which contain one term for the diffusion of CO 2 /O 2 into/out of leaves through stomata, and one term for the depletion of CO 2 and production of O 2 because of photosynthesis and one term for the production of CO 2 and depletion of O 2 due to respiration. Assuming that diffusion of these gases is governed by Fick's law (Fick, 1855), for the diffusion term we get the fluxes Here g w is the stomatal conductance to water in mol m 2 s , C A is the atmospheric CO 2 concentration in mmol mol and O A is the atmospheric O 2 concentration in mmol mol . The ratios of the diffusivities in air of water, carbon dioxide and oxygen mean that g w 1:6 is the stomatal conductance to carbon dioxide and g w 1:25 is the stomatal conductance to oxygen (Haynes et al., 2016). The fluxes J c and J o both have the unit mol m 2 s . Dividing them by the thickness of the leaves d L in m we get the rate of change of gas molecules per unit volume, which has units mmol m 3 s . Since we denote concentrations in mmol mol , we want the rate of change per mol, instead of per unit volume. According to the ideal gas law, n = PV RT ,where V is a volume of gas in m 3 , P the pressure of that gas in Pa, T the temperature of that gas in K and R the ideal gas constant in J Kmol So by multiplying the fluxes by RT a PD L we get the rate of gas concentration change in our preferred units, where T a is the atmospheric temperature in K. Including terms for the photosynthetic assimilation rate A and the dark respiration rate R d , both in mmol m 2 s , we get Here A is given by equation 2. The gas concentrations in leaves typically reach equilibrium conditions in less than a minute. Since typically OpenSimRoot timesteps are about a tenth of a day, we will make a quasi steady-state assumption. At steady-state, dC m dt = dO m dt = 0 and equations 12 and 13 give In C 4 photosynthesis, the mesophyll exchanges gases with the atmosphere and CO 2 is then transported to the bundle sheath. Therefore, we model mesophyll and bundle sheath gas concentrations as connected reservoirs, using the differential equations: where C s and O s are the bundle sheath CO 2 and O 2 concentrations in mmol mol , R m and R s are the respiration rates in mesophyll and bundle sheath cells in mmol m 2 s , respectively, g s is the conductance of the mesophyll-bundle sheath interface to CO 2 mol m 2 s , 0.047g s is the mesophyll-bundle sheath interface conductance to O 2 (Von Caemmerer, 2000), V P is the PEP carboxylation rate in mmol m 2 s , and A r is the assimilation rate, not including dark respiration, in mmol m 2 s . The equations concerning mesophyll concentrations have terms for gas exchange with both the atmosphere and the bundle sheath cells, a respiration term and a term for the PEP carboxylation, while the equations concerning bundle sheath concentrations have a term for the gas exchange with the mesophyll cells, a term for respiration and a term for photosynthesis. PEP carboxylation rates are given by Here V pmax in mmol m 2 s is the maximum PEP carboxylation rate, K P in mmol mol the Michaelis constant for PEP carboxylation and V pr in mmol m 2 s the PEP regeneration rate. As for C 3 photosynthesis, gas concentrations reach equilibrium values very quickly as compared to the timescale relevant to OpenSimRoot and relative to changes in environmental conditions. Thus, we make a quasi steady-state assumption and we get:

Stomatal conductance
For the stomatal conductance to water, g w , we use a slightly modified version of the Ball-Berry-Leuning model, described in (Leuning, 1995), which states that Here m is an empirical constant, S w (y c ) is the water stress factor defined in Section 2.1, Г in mmol mol is the CO 2 compensation point with dark respiration, VPD in Pa the vapor pressure deficit, VPD r in Pa a reference vapor pressure deficit and g w0 in mol m 2 s the residual conductance. For C 3 photosynthesis, Г is given by From (Von Caemmerer, 2000) we get the C 4 expression for Г, which is Here Г s in mmol mol is the bundle sheath CO 2 concentration at the compensation point, which is given by

Leaf temperature
Like many other aspects of plant development and function, photosynthetic assimilation rates are temperature dependent (Duncan and Hesketh, 1968;Lin et al., 2012). In order to make accurate predictions in the context of climate change, it is important that we can model the effects of elevated temperature and increasing atmospheric CO 2 concentrations. This requires us to model the temperature dependence of the kinetic parameters that our models rely on (Bowes, 1991;Long, 1991;McMurtrie and Wang, 1993;Walcroft et al., 1997). For most parameters, we follow (Von Caemmerer, 2000) and use an Arrhenius function of the form where T L is the leaf temperature in K, T ref the reference temperature, 298.15 K (25°C), P(T L ) is the value of the parameter at T L , P 25 is the parameter value at the reference temperature of 298.15 K (25°C), E in J mol is the activation energy and R is the universal gas constant. We model the temperature 2S c=o and V cmax with equation 29. For other parameters, the peaked Arrhenius function of the form is used because it better represents the temperature dependence of that parameter. Here D in J mol is the deactivation energy and S in J Kmol is called the entropy factor. We use equation 30 to model the temperature dependence of J max , V pmax and g s . To calculate leaf temperature T L , we use an adaptation of the leaf energy balance as described in (Nobel, 2009). The terms that are described there as contributing to the leaf energy balance are: absorption of solar irradiation A S ; absorption of infrared radiation from the surroundings A IR ; emission of infrared radiation e IR ; heat conduction and convection H c ; heat loss due to evapotranspiration H t ; photosynthesis; metabolic processes. The contribution of photosynthesis and metabolic processes to the leaf energy balance is typically on the order of a few Watts, while other terms such as the solar irradiation are several hundreds of Watts, so to simplify the model we omit the photosynthesis and metabolic processes from the energy balance. Due to the low specific mass of leaves (approximately 0:2 kg m 2 ), their specific heat is low and a net energy balance of a few Watts is enough to increase the temperature by one kelvin in less than 5 minutes. Because of this, we assume the leaf is in steady state and the energy balance E is zero: The absorption of solar irradiation per unit area A s in W m 2 is equal to where a is the absorptance of the leaf over the whole solar spectrum, r is the fraction of solar irradiation reflected by the surroundings and I in W m 2 is the solar irradiation. Assuming there are nonearby objects that reflect a lot of sunlight towards the field in which we are modelling crops and using a symmetry argument (for each photon reflected from another leaf, we can assume the leaf reflects a photon itself), we set r = 0. So that Using the Stefan-Boltzmann law, the absorption of infrared radiation per unit area A IR in W m 2 is Here a IR is the infrared radiation absorptance of the leaves, s in W m 2 K 4 is the Stefan-Boltzmann constant, T surr in K is the temperature of the surroundings and T sky in K is the effective temperature of the sky. Here we assume that the underside of the leaves absorb infrared radiation coming from the surroundings while the upper side absorbs infrared radiation coming from the sky. Note that T sky is not the actual temperature of the sky, but the temperature a blackbody that emitted as much radiation as the sky would have. We assume that T sky = T a − 40: Here T a is the air temperature in K. A blackbody with absorptance a will also have emissivity a so the emission of infrared radiation by the leaves, e IR in W m 2 , is equal to The factor 2 in this equation 37 comes from the fact that the leaves emit infrared radiation from both sides. The heat loss due to conduction into the air boundary layer around leaves and then convection away from the leaves, H c in W m 2 , is equal to Here K air (T a ) in W mK is the thermal conductivity of the air, d bl in m the thickness of the boundary layer, d in m the characteristic length of the leaf, v (T a ) in m 2 s the kinematic viscosity of the air and v in m s the wind speed. The factor 2 in equation 38 is because heat is conducted away from both sides of the leaf. For the range of air temperatures we are concerned with, 273.15 K to 323.15 K (0 to 50°C), we use, from (Nobel, 2009), the following approximations K air T a ð Þ = 0:0243 + 0:00007 T a − 273:15 ð Þ , n T a ð Þ = 1:415 + 0:09 T a − 273:15 ð Þ ð Þ · 10 −5 : The heat loss due to evapotranspiration, H t in W m 2 is equal to where J v in mol m 2 s is the transpiration rate and H vap (T a ) in J mol is the heat of vaporization of water, which is equal to Putting this together and setting the leaf energy balance equal to zero we obtain This is a quartic equation in T L , the leaf temperature, which we solve numerically using the Newton-Raphson method. As our initial guess we take T a , the atmospheric temperature, which should be somewhat close to T L .

Scaling from leaf to canopy
The photosynthesis, leaf gas concentration, stomatal conductance and leaf temperature models described in previous sections are all models for individual leaves. However, in a canopy, not every leaf is subjected to the same conditions. Some leaves are in direct sunlight, while others will be partially or completely shaded by other leaves. A number of different models have been proposed to scale up from the leaf scale to the canopy scale. The simplest, big leaf models, approximates the entire canopy as a single leaf (Amthor, 1994;Lloyd et al., 1995). These models tend to overestimate photosynthesis rates and require empirical extinction factors to make them produce accurate results. More sophisticated models use multi-layer approaches, which divide the canopy up in layers, applying leaf-scale models to each separately (de Wit, 1965;Duncan et al., 1967;Waggoner and Reifsnyder, 1968;Lemon et al., 1971;Chang et al., 2018). We use a third approach, the sun/shade model described in (De Pury and Farquhar, 1997), which divides the canopy up into sunlit and shaded leaf area fractions, producing results with similar accuracy to the multi-layer model (De Pury and Farquhar, 1997) while also being relatively simple. This is the model we chose to represent our canopy. In this model the sunlit and shaded leaf area fractions are calculated each time step based on leaf area index and the angle of the sun with respect to the canopy. We integrate this with the previously described models by calculating stomatal conductance, leaf gas exchange, leaf temperature and photosynthesisrates separately for the shaded and sunlit leaf areas.
In order to accurately simulate the incoming solar radiation required by the multi-layer model described in (De Pury and Farquhar, 1997), the model described in (Michalsky, 1988) was implemented. This allows OpenSimRoot to simulate the amount of incoming solar radiation to high precision. It also allows us to calculate related quantities like the day length.

Nitrogen limitations on photosynthesis
Leaf nitrogen content is an important factor determining photosynthetic assimilation rates (Sinclair and Horie, 1989). Using existing OpenSimRoot models for minimum and optimum nutrient concentrations, which have to be specified for every plant organ, and the total plant nitrogen content, we estimate leaf nitrogen content leaf N : where leaf min and leaf opt are the minimal and optimal leaf nitrogen contents, plant min and plant opt are the minimal and optimal plant nitrogen contents and plant N is the total plant nitrogen content. We assume the electron transport rate J max and carboxylation rate V max are limited by leaf nitrogen content through the folowing equations: where n 1 and n 2 are proportionality constants (Kull and Kruijt, 1998) and N leaf is the leaf nitrogen concentration (not content as in equation 45). By substituting J max into equation 5 and V cmax in equation 3 (or equation 7 in the case of C4 photosynthesis), nitrogen limitations are taken into account.

Carbon reserves
OpenSimRoot contains a relatively simple model for the accumulation and degradation of starch carbon reserves. From the available carbon, C a in g day , from photosynthesis and the seed, carbon costs such as respiration and root exudation, C c in g day , are first deducted. Potential growth rates of different plant organs together determine the carbon sink for growth, C g in g day . Any remaining carbon after allocating to growth is added to the carbon reserves (a pool representing the carbon in starch and other metabolites). The allocation rate of carbon to reserves is denoted by C r in g day (note that if carbon is remobilized from the reserves, C r <0), while the current amount of carbon in the reserves is denoted by C p in g. So in most cases If not enough carbon is available from photosynthesis and the seed for carbon costs and growth, so if C a < C c +C g , carbon from the reserves is remobilised. The plant cannot remobilize more than half of the current reserves per day, so C r > − C p 2day . If carbon from photosynthesis, the seed and reserves together is not enough to match the carbon costs and growth, so if C a -C r < C c +C g , the allocation of carbon to growth will be less than the amount required for potential growth. The actual amount of carbon available for growth, C h , will then be C h = C a -C r -C c < C g . This will lead to a reduction in growth rates, thus controlling growth through carbon availability.
This carbon reserve allocation model works well with a photosynthesis model based on light-use efficiency which averages out photosynthesis over the full 24 hours in a day. However, since we are now adding models for simulating the day-night cycle, we need to allow for carbon to be stored during the day and then released during the night. The model has to satisfy these important requirements: • Carbon remobilization from the reserves at night should be high enough to allow for uninterrupted growth if water and nutrients are plentiful. • Allocation of carbon to the reserves during the day should be high enough such that there is enough in reserves to match the non-growth carbon costs during the night. • If photosynthesis rates over a 24-hour period are lower than the carbon needed for maintenance and growth then growth rates should be reduced.
Based on our requirements, we implement the following modifications to the model of carbon reserves.
• At night, the maximum rate of starch remobilization is determined at the start of the night and set such that at most 95% of the reserves are remobilized during the night. So C r > −0:95 C p t n where t n in d is the duration of the night.
• During the day, if carbon reserves at the end of the day are expected to be below 110% of the total carbon needed during the night, growth is reduced. More precisely, if C p + C r t s ≤1.1C c N, where t s in d is the time until sunset, then C r is bounded from below by min (C a − C c , 1:1C c N−C p t s ).
These two modifications are simplifying model assumptions based on an analysis of Arabidopsis carbon reserves under photoperiods of different length which showed that starch reserves were remobilized at fairly constant rates which depended on the length of the photoperiod (Sulpice et al., 2014). It was also found that starch accumulation was more rapid in shorter photoperiods. While the diurnal dynamics of growth differ between monocots and dicots, these differences are not related to altered carbohydrate availability for growth (Poire et al., 2010). Since OpenSimRoot does not simulate inherent (i.e., not resource limited) daily variation in growth rate, nor does it simulate temperature-dependent growth rates, there is no need to make a distinction between monocots and dicots in the carbon reserve allocation model.

OpenSimRoot root finder
At any given time, in order to calculate the photosynthetic assimilation rates, we have to find a solution to the system of equations 2, 3, 4, 14, 15 and 25 (in the case of C4 photosynthesis this is 21, 22, 23, 24, 20, 7, 8, 9 and 25). It has been shown that when solving the coupled photosynthesis and stomatal conductance equations using a fixed point method there is no convergence in many cases (Sun et al., 2012). In order to guarantee convergence and arrive at the correct solution in that case, it has been suggested to use the Newton-Raphson method (Sun et al., 2012). Ideally, we would use a multivariate Newton-Raphson method to solve the above system of equations. However, this is only possible within the OpenSimRoot engine if the whole system of equations and relevant solutions are solved within a single submodel. Doing this would make it much more complicated for future users to modify the equations, choose different models for one of the relevant state variables or add new models coupled to the equations above, which is what the modular structure of the OpenSimRoot engine was designed to facilitate. By keeping models separate, we allow users to pick and choose between different models for variables like photosynthesis, leaf gas exchange, leaf temperature or stomatal conductance. Because of this, we have implemented a root finder based on the Newton-Raphson method which works within the constraints of the OpenSimRoot engine. Suppose we have a system of equations x n = F n x 1 , :::, x n ð Þ , where x 1 could for example be the stomatal conductance and F 1 equation 25, the expression for stomatal conductance. We start from initial values x 1,0 , …, x n,0 . The root finder will then use the Newton-Raphson method to find an x n,* for which F n x 1,0 , :::, x n−1,0 , x n, * − x n, * < Є n : Here Є n is a tolerance. Then, we want to find a x n-1,* which satisfies F n−1 x 1,0 , :::, x n−1, * , x n, * − x n−1, * < Є n−1 : We can again do this using the Newton-Raphson method but at each step will update x n,* . More concretely, the root finder does one step of the Newton-Raphson method to find x n-1,1 and finds a new x n,* that now satisfies F n x 1,0 , :::, x n−1,1 , x n, * − x n, * < Є n : After a number of steps we find an x n-1,* which satisfies F n x 1,0 , :::, x n−2,0 , x n−1, * , x n, * − x n−1, * < Є n−1 : Now we use the same process to find x n-2,* which satisfies F n x 1,0 , :::, x n−2, * , x n−1, * , x n, * − x n−2, * < Є n−2 , again updating x n-1,* and x n,* at each iteration step. Following this iterative procedure for all variables we eventually find x 1,* ,x 2, * ,…,x n,* which satisfy F 1 x 1, * , :::, x n, * − x 1, * < Є 1 (56) ::: F n x 1, * , :::, x n, * − x n, * < Є n : This is the simultaneous solution to the system of equations making up the photosynthesis model. If there is no convergence in 40 steps, the value from the previous timestep will be used (in case this happens on the first timestep, a predetermined default value will be used).

Simulation setup
In this paper, we test the utility of several variants of the "Steep, Cheap, and Deep" (SCD) ideotype against a reference root phenotype of maize under rainfed and drought conditions in three distinct maize-producing pedoclimatic environments. The SCD phenotypes are different from the reference phenotype in 4 phene states: • The number of axial roots from different classes (i.e., seminal, nodal, and brace) is reduced from 89 to 45. • The angles of all axial root classes (i.e., seminal, nodal, and brace) are 20 degrees steeper compared to the reference phenotype. • The lateral root branching density (LRBD) is reduced by 50%. • The maximum root cortical aerenchyma (RCA) volume is increased from 20% to 40%.
In addition to simulating the reference and SCD phenotype, we simulated every combination of the reference and SCD values for the above 4 phenes resulting in 16 different combinations. We will refer to these different phenotypes as follows: • All phenes have the reference value: Reference phenotype. • All phenes have the reference value except one which has the SCD value: The phene change from the reference value (so "Fewer axial roots", "steeper axial roots", "reduced lateral root branching density", "higher root cortical aerenchyma"). • Two phenes with the reference value, two with the SCD value: The phene values which are not the reference value (so for example: "Fewer, steeper axial roots"). • One phene with the reference value, three with the SCD value: SCD + the phene value which has the reference value (so for example: "SCD + more axial roots" if the first phene from the list above has the reference value). • All four phenes have the SCD value: Steep, cheap, deep (SCD).
We test the utility of these 16 root phenotypes in 6 different pedoclimatic environments including rainfed and vegetative terminal drought conditions. It should be note that OpenSimRoot is a mechanistic model which is not calibrated to match specific experimental data. As a heuristic model, it is designed to test the adequacy of a hypothesis, provide explanations and identify knowledge gaps. Like any model, it excludes many factors that may affect the actual performance of plants in the field, however it is never calibrated to force agreement with empirical results, instead using empirically measured root phenotypes and models to generate predictions. Simulation studies with OpenSimRoot aim to uncover qualitative relationships and quantifying the effects of root system phenotypes, rather than providing exact estimates of plant performance.
In every simulation in this study, a single maize plant is simulated representing an individual within a monoculture stand with a between-row spacing of 50, within-row spacing of 25, and soil depth profile of 150 cm, corresponding to a planting density of 8 plants per m 2 and sowing depth of 5 cm. Roots from neighboring plants are simulated by mirroring the roots at the boundary back into the column. To accurately mimic the rainfed maize agrosystem, actual soil and climate datasets for three different maize growing locations, namely Boone County (Iowa, USA -42°49.4'N 93°43'45.9"W), Zaria (Kaduna, Nigeria -11°1'N 7°37'E), and Tepatitlań (Jalisco, Mexico -20°5 2'N 102°43'W) are used. For simplicity, these sites are referred to as Iowa, Zaria, and Jalisco. The soils in these sites have been identified as Mollisol -Udoll (Iowa), Alfisol -Ustalf (Zaria), and Andosol -Vitrand (Jalisco). Soil profile data, including soil texture, soil bulk density, and relative water content at -33 and -1500kPa, from the ISRIC soil hub database (https://data. isric.org/) is used for Iowa and Jalisco, and for Zaria a published dataset is used (Beah et al., 2020). These soil profile datasets are implemented to estimate the corresponding van Genuchten parameters using the online freeware Rosetta3 (https://www. handbook60.org/rosetta/). The weather data for the year 2019 over the main growing season for all three locations come from the POWER LaRC database (https://power.larc.nasa.gov/dataaccess-viewer/). For the drought scenarios, we simulate a terminal vegetative stage drought by setting precipitation rates to zero after a certain number of weeks. This time point is chosen such that the reference phenotype would see a reduction in shoot dry weight of approximately 50% as compared to the rainfed scenario. For Iowa, the terminal vegetative drought starts on day 21, while for the Zaria and Jalisco it starts on day 28.
The OpenSimRoot code is available at https://gitlab.com/ rootmodels/OpenSimRoot. Instructions on how to acquire the correct executable and reproduce the simulations done in this study can be found at the OpenSimRootpaper files repository. See (Schäfer et al., 2022b) for a comprehensive introduction on running OpenSimRoot. Supplementary File 1 contains an input file with all parameters corresponding to the new capabilities. Supplementary File 2 contains a summary of the simulated results.

Estimating phene interactions
The actual (predicted) shoot biomass in response to the drought of each integrated SCD phenotype is compared to their corresponding expected shoot biomass response to quantify interactions among root phenes. The expected shoot biomass for each SCD phenotype is derived using the additive null model (Coteé t al., 2016) and the shoot biomass response of the reference maize root phenotype. The shoot biomass responses of the individual phene states are calculated by subtracting the reference shoot biomass from each SCD phenotype that differs in just one phene state. The resulting shoot biomass responses of individual phene states are then added to the reference shoot biomass to get the expected additive shoot biomass corresponding to a combination of phene states. The variance of the expected additive shoot biomass is calculated by summing the variances of the relevant quantities (the variance of a sum of distributions is the sum of variance). From this variance we calculate the standard error by taking 5 as the population size, matching the number of repetition simulations we run for each set of inputs. The actual response greater than, equal to, and less than the expected response respectively corresponds to the synergistic, additive, and antagonistic relationship between the root phenes. More details on this can be found in (Ma et al., 2001;Ajmera et al., 2022).

Pedoclimatic conditions alter the utility of root phenotypes
The results highlighted in this study demonstrate the new capabilities that mechanistically capture processes underlying shoot physiology in OpenSimRoot_v2 (Figure 1). The architecture of maize root phenotypes ( Figure 2) and in turn, their impact on shoot dry weight varied in response to different pedoclimatic conditions, including drought (Figure 3). In the rainfed Iowa environment, the reference phenotype and the high aerenchyma phenotype had the greatest shoot dry weight with just under 25 g. In response to drought in Iowa, the shoot dry weight of the reference phenotype was reduced to 13 g. However, phenotypes with fewer, shallow axial roots and reduced lateral root branching density had the greatest shoot dry weights under drought in Iowa (approximately 20 g). Steep axial roots were not associated with high shoot dry weight in Iowa drought conditions and the amount of RCA had little effect on shoot dry weight in these conditions.
In Zaria, the reference phenotype had intermediate performance in both the rainfed and drought conditions, with 19.75 and 9.5 g of shoot dry weight. In rainfed conditions, two phenotypes with shallow axial roots and abundant aerenchyma had shootdry weights smaller than 10 g, indicating that rainfed Zaria conditions were more stressful than the rainfed Iowa conditions. In Zaria, reduced lateral branching density was associated with increased shoot dry weight in combination with most other phenesunder drought conditions. In contrast to Iowa, reducing the number of axial roots without any other changes in root phenotypes was slightly detrimental in terms of shoot biomass gain in response to drought in Zaria. In Jalisco, both the reference and SCD phenotypes were near the greatest shoot dry weight in both rainfed and drought conditions. The phenotypes with reduced lateral branching density (LRBD), few axials plus reduced LRBD, reduced LRBD plus abundant aerenchyma, and shallow axial roots were the best performers in both environments. Except for rainfed Iowa, across all the environments the phenotypes with steep axial roots plus either level of aerenchyma had the smallest shoot dry weight.

'Steep, Cheap, Deep' phenotype reduces the carbon cost of soil exploration
Under rainfed conditions, the SCD phenotype had notably lower root carbon costs (structural carbon, respiration, and exudation) than the reference phenotype, with a reduction of approximately 30% in Zaria (Figure 4). Under drought, the difference between root carbon costs of the reference and SCD phenotypes decreased in the Iowa and Zaria environments and stayed the same in the Jalisco environment. Decreasing the axial root number reduced root carbon costs the most in all environments except the Jalisco drought environment. The effect of other single phene changes depend more on the environment. Phenotypes with few axial roots combined with abundant aerenchyma and large lateral root branching density had small root carbon costs in most environments.
The carbon use efficiency (water acquired per unit of carbon invested in roots) of the SCD phenotype was greater than that of the reference phenotype in all environments, (although the difference was very small in Jalisco ( Figure 5). From the individual changes, fewer axial roots, reduced lateral root branching density, and greater aerenchyma formation improve carbon use efficiency, while steeper axial roots reduce carbon use efficiency. In Iowa, the carbon use efficiency was greater than in the other two environments by about 0.5 liters per gram of carbon.

'Steep, Cheap, Deep' phenotype facilitates rooting in deeper soil horizons
The SCD phenotype, with either aerenchyma level, has the largest fraction of root length deeper than 50 cm ( Figure 6). In Iowa, 1% or less of roots in the reference phenotype were deeper than 50 cm, versus more than 8% for the SCD phenotype. In Zaria, 10% of the root length was below 50 cm for the reference phenotype, which increased to 20% in the rainfed environment and 24% in the drought environment for the SCD phenotype. In Jalisco, about 17% of root length was below 50 cm for the reference phenotype, which increased to ca. 30% for the SCD phenotype. Reducing axial root number lead to a bigger increase in deep roots than making the axial root angle steeper, as does reducing lateral root branching density in the majority of cases. Aerenchyma formation had little effect on the fraction of deep roots.

'Steep, Cheap, Deep' phenotype encounters lower penetration resistance in diverse pedoclimatic conditions
There were major differences in the depth profiles of soil penetration resistance encountered by roots among the 6 environments ( Figure 7). In general, rainfed environments had less soil penetration resistance than drought environments: the difference between them was often several thousand kPa. In rainfed Iowa soil, penetration resistance sharply increased between 20 and 45 cm. Under drought, Iowa soil had the Visualization of reference (panels A, C, E) and SCD (panels B, D, F) maize root phenotypes in the Iowa (panels A, B), Jalisco (panels C, D) and Zaria (panels E, F) environments without drought stress at 42 DAG.
greatest penetration resistance of all the soils, peaking at more than 14000 kPa at 45 cm. The soil penetration resistance of the rainfed Zaria soil increased from around 900 kPa at the soil surface to 5000 kPa at 150 cm. In the Zaria drought environment, the soil penetration resistance was greater at every depth and reached its maximum values between 50 and 80 cm deep, after which it declined again. In the Jalisco rainfed environment, the reference phenotype encountered greater soil penetration resistance than the SCD phenotype, the soil penetration resistance encountered by the reference phenotype is similar to what both phenotypes encounter in the Jalisco drought environment.

'Steep, Cheap, Deep' phenotype retains more water at the root surface
Plants with the SCD phenotype had more water near their roots at most depths than plants with the reference phenotype ( Figure 8). In the rainfed Iowa and Zaria environments, there was more water near the roots in the topsoil than there was in deeper soil layers. In drought environments, there was less water near the roots than in rainfed environments and there was relatively little water near the roots in the topsoil and more in deeper soil strata. Shoot dry weight of all the maize phenotypes considered after 42 days in the 6 different environments. The error bars indicate standard deviation, n (replicates) = 5. LRBD means lateral root branching density, RCA means root cortical aerenchyma and SCD refers to the steep, cheap, deep phenotype.

'Steep, Cheap, Deep' phenotype enhances carbon deposition in deeper soil horizons
The SCD phenotype had notably greater carbon deposition (carbon in the form of root biomass or exudates) below 50 cm in the soil than the reference phenotype in all 6 environments ( Figure 9). In Iowa, the SCD phenotype deposited approximately 8 times more carbon deep in the soil than the reference phenotype under both rainfed and drought conditions. In the Zaria soil, the SCD phenotype deposited 1.6 times more deep carbon in the rainfed case and 2.2 times more in the drought case, while in the Jalisco soil, the difference was a factor 2 in both cases. Besides the SCD phenotype, a few other phenotypes had similar or even greater carbon deposition. The performance of the SCD phenotypes with less aerenchyma formation or shallow axials were similar in most environments. This was similar to the performance of the phenotypes with few steep axial roots and those with low lateral root branching density in several environments. Except in the rainfed Iowa environment, there was at least one SCD variant phenotype that had both greater shoot dry weight and greater carbon deposition below 50 cm than the reference phenotype ( Figure 10). In the rainfed Iowa environment, the reference phenotype had the greatest shoot dry weight and least deep soil carbon deposition while the SCD phenotype had one of the Total carbon invested in roots (biomass, exudation and respiration) of all the maize phenotypes considered after 42 days in the 6 different environments. The error bars indicate standard deviation, n (replicates) = 5. LRBD means lateral root branching density, RCA means root cortical aerenchyma and SCD refers to the steep, cheap, deep phenotype. greatest deep soil carbon depositions and smallest shoot dry weights. In almost all environments, steeper axial roots by themselves notably reduced shoot dry weight for a relatively small increase in deep soil carbon deposition.

'Steep, Cheap, Deep' phenotype improves deep soil water capture
Plants with the SCD phenotype had less water capture in the upper 40 cm of soil and greater water capture in deeper soil strata ( Figure 11). Plants with the SCD phenotype also acquired water from deeper in the soil. The difference in water capture from the subsoil was greater in drought environments than in rainfed environments. In some environments, negative total water uptake in the topsoil highlights the occurrence of hydraulic lift.

Phene synergism influences plant performance and varies with pedoclimatic conditions
We observed substantial synergism (i.e., greater benefit than expected from simply additive interaction) among root phene states for plant performance in all environments except rainfed conditions in Iowa (Figure 12). The magnitude of interactive effects was Carbon use efficiency (water uptake per unit of carbon invested in roots) of all the maize phenotypes considered after 42 days in the 6 different environments. The error bars indicate standard deviation, n (replicates) = 5. LRBD means lateral root branching density, RCA means root cortical aerenchyma and SCD refers to the steep, cheap, deep phenotype. substantial, ranging from antagonism of -36% (SCD plus shallow axial roots in rainfed Zaria, middle left panel, ACD in Figure 12) to synergism of +226% (SCD plus greater lateral root branching density in Jalisco under drought, bottom right panel, ABD in Figure 12), averaging +30% over the additive effects. In Zaria and Jalisco, 6 integrated phenotypes showed strong synergism in multiple environments, of which all 6 include steeper axial root growth angles, 4 include fewer axial roots and 4 include reduced lateral root branching density. In Iowa under drought, two integrated phenotypes displayed synergism, both containing the combination of steeper axial root growth angles and increased aerenchyma formation. Synergism was greater under drought than under rainfed environments, averaging +39% under drought compared with +24% under rainfed conditions in Zaria, and +74% under drought compared with +42% under rainfed conditions in Jalisco, all relative to the additive interactions ( Figure 12). Synergism was smaller and less common in the Iowa environment.

Discussion
The implementation of established biophysical-biochemical models to mechanistically capture various shoot processes in the Fraction of the total root length below 50 cm deep in the soil of all the maize phenotypes considered after 42 days in the 6 different environments. The error bars indicate standard deviation, n (replicates) = 5. LRBD means lateral root branching density, RCA means root cortical aerenchyma and SCD refers to the steep, cheap, deep phenotype.

FIGURE 8
Volumetric water content near the roots on day 42 for 12 individual simulations, averaged over 1 cm layers. Soil penetration resistance felt by the roots on day 42 for 12 individual simulations, averaged over 1 cm layers. Schäfer et al. 10.3389/fpls.2022.1010165 Frontiers in Plant Science frontiersin.org functional-structural plant/soil model OpenSimRoot has led to the development of one of the most feature-rich root/soil modeling platforms, OpenSimRoot_v2, available in the literature to date. Using the OpenSimRoot_v2 model, we evaluated the utility of SCD root phenotypes under rainfed and drought conditions over three distinct pedoclimatic environments to evaluate the relationship of root phenotypes with maize growth and soil carbon deposition under drought. Our results show that parsimonious root phenotypes with fewer axial and lateral roots improve shoot dry weight notably in terminal drought conditions, as hypothesized (Passioura, 1983;Lynch, 2013). We observe that the SCD phenotype spends less carbon on the root system without negatively affecting water capture. This results in greater carbon use efficiency (water acquired per unit of carbon invested in roots) that in turn leads to greater shoot biomass. The increase in root length at depth allows plants with the SCD phenotype to access deeper water, and distribute water capture over more soil strata, so roots encounter softer, wetter soil, and increase deep soil carbon deposition. Reduced axial root number has been suggested as a target for breeding more drought-resistant plants (Lynch, 2018). Maize plants with fewer crown roots had greater shoot biomass under drought in mesocosms and the field (Gao and Lynch, 2016). Reducing axial root number increased shoot dry weight in the Iowa drought environment but it slightly reduced or had little effect on shoot dry weight in the other environments ( Figure 3). Reducing axial root number decreased root carbon costs and increased root length and Carbon depositions below 50 cm deep in the soil of all the maize phenotypes considered after 42 days in the 6 different environments. The error bars indicate standard deviation, n (replicates) = 5. LRBD means lateral root branching density, RCA means root cortical aerenchyma and SCD refers to the steep, cheap, deep phenotype. relative root length below 50 cm (Figures 4, 6). This reduced root carbon cost and a deeper root system were not associated with an increase in shoot dry weight.
Having steeper axial roots by itself was detrimental in all 6 of the simulated environments ( Figure 2). This result highlights the importance of capturing water available in the topsoil, even in drought environments where more water is found in the subsoil. Shallow root systems are subject to greater interplant competition for topsoil resources like phosphorus (Lynch and Brown, 2001;Rubio et al., 2001) while steeper root systems experience greater intraplant competition for mobile resources like nitrogen and water (Nord et al., 2011;Trachsel et al., 2013;Ajmera et al., 2022). Steeper axial root angles are associated with deeper root systems (Oyanagi et al., 1993;Uga et al., 2015), but the benefit of steeper axial roots depends on environmental factors such as the availability of subsoil water and the amount of Carbon depositions below 50 cm deep in the soil versus shoot dry weight of all the maize phenotypes considered after 42 days in the 6 different environments. LRBD means lateral root branching density, RCA means root cortical aerenchyma and SCD refers to the steep, cheap, deep phenotype.
in-season precipitation (Manschadi et al., 2008). Root systems with both greater root length density in intermediate and deep soil layers and deeper maximum depth are candidates for increased subsoil water uptake and thus drought resistance (Hund et al., 2009;Wasson et al., 2012). Making the axial root angles of the reference phenotype steeper increased root length below 50 cm, but markedly less than doing this in combination with reducing axial root number and lateral branching density, which highlights the importance of interactions between root phenotypes and the environment in determining drought resistance (Ludlow and Muchow, 1990;Richards et al., 2002). While increased root length in the subsoil was beneficial in the terminal drought scenario, large and shallow root systems were beneficial in some simulated environments. Capturing water in the topsoil, which is of course where precipitation enters the soil, will be even more important for crops that have to compete with weeds and intercrops or cover crops for limited water resources.
Reducing lateral root branching density increased shoot dry weight in all environments except the rainfed Iowa environment (Figure 3). Since competition among neighboring lateral roots for water is significant, reducing lateral branching density reduces root system carbon costs much more than it reduces water capture. In fact, the reduction in lateral root carbon costs due to reduced lateral branching density increased photosynthate available for axial roots, which were then able to grow deeper, increasing root depth ( Figure 6). These results are consistent with the report that reduced lateral root branching density was associated with deeper rooting, better capture of subsoil water, and better plant water status, growth and yield in maize under drought stress in two field environments and in greenhouse mesocosms . Similarly, maize lines with reduced lateral root branching density had deeper rooting, greater plant N content, photosynthesis, growth, and yield under low N stress in two field environments and in greenhouse mesocosms .
Increasing aerenchyma formation had little effect on shoot dry weight in all environments (Figure 3). Likewise, for most intermediate phenotypes to which increased aerenchyma formation was added, shoot dry weight did not change much. This contrasts with field studies where increased aerenchyma formation was associated with substantially greater maize yield under drought in the field (Zhu et al., 2010;Chimungu et al., 2015). This discrepancy is perhaps due to the narrow range of increase in aerenchyma formation (i.e., from 20% to 40%) considered in this study compared to the field experiments wherein aerenchyma formation was observed to vary from 0 to 50% in maize under drought with significant stress and environmental plasticity (Schneider et al., 2020). Aerenchyma influences plant performance by reducing root respiration and remobilizing tissue nutrients, which are explicitly simulated in the model. Since nutrient constraints on growth are not included in this drought study, the benefit of aerenchyma formation is less than that observed in response to suboptimal nutrient availability (Postma and Lynch, 2011).

FIGURE 11
Relative root water uptake after 42 days for 12 individual simulations, averaged over 3 cm layers. There is negative uptake in the upper soil layers because of hydraulic lift.
Our observation that phene synergism is an important component of the fitness of integrated root phenotypes supports previous reports of substantial phene synergism for soil resource capture in silico and in planta. For example, increased root hair length, root hair density, root hair proximity to the root tip, and number of trichoblast files increased P acquisition in Arabidopsis roots by 371% more than their additive effects in silico (Ma et al., 2001). In simulated maize roots, Root cortical aerenchyma increased P acquisition 2.9 times more in phenotypes with increased lateral branching density (Postma and Lynch, 2011).
Multiple synergies are evident in the interactions of axial root phenotypes of common bean for the acquisition of N and P in silico (Rangarajan et al., 2018). In planta, metaxylem anatomy and root depth are synergistic for drought adaptation in contrasting Phaseolus species (Strock et al., 2021), and shallow root growth angles and long root hairs aresynergistic for P capture in common bean (Miguel et al., 2015). Here, we attribute synergism between fewer axial roots and steeper axial root growth angles to reduced intraplant competition for water capture. Steeper root growth angles increase subsoil exploration but also position root axes Shoot dry weight of phenotypes as expected by adding the changes resulting from altering individual phenes compared to the actual shoot dry weights of these phenotypes. The black line indicates the shoot dry weight of the reference phenotype. The error bars indicate standard error. Cases where the expected and actual shoot dry weight are significantly different are marked by an asterisk on the x axis label. closer together, thereby increasing competition among root axes for soil resources, especially mobile resources such as water (Dathe et al., 2016). Two synergistic integrated phenotypes combine steeper axial root growth angles with reduced lateral root branching density, which was synergistic in 4 of the 6 environments. Reduced lateral root branching, as with reduced axial root number, would reduce competition among roots for soil water. In this context we note that slight water stress occurs under our rainfed Zaria and Jalisco scenarios (Figure 3), meaning that root phenotypes that improve water capture should improve plant growth under both rainfed and drought environments. We interpret the smaller degree of synergism observed in the Iowa environments to the notably greater bulk density of the Iowa soil, which inhibits soil penetration by small diameter roots, reducing intraplant competition for water capture.
While shoot dry weight under drought correlates very strongly with water uptake, we also found a strong relationship between carbon use efficiency (water uptake per unit of carbon invested in roots) and shoot dry weight. A similar metric, root system efficiency (transpiration per unit leaf area per unit of root mass) was proposed as breeding target for drought resistant maize (van Oosterom et al., 2016). By reducing the amount of carbon used for growing and maintaining the root system, a plant frees up carbon which can be used to grow a bigger shoot or for meeting metabolic needs during drought stress (Passioura, 1983;Bolaños et al., 1993;Lynch, 2018). If the carbon required to grow extra roots do not result in enough additional water uptake to produce more carbon, the plant is effectively moving biomass underground, reducing its yield potential. Three of the 4 phene states included in the SCD phenotype reduce root system carbon costs: Fewer axial roots, reduced lateral branching density and greater root cortical aerenchyma formation. Noting that reducing the axial root number and lateral branching density increase carbon use efficiency more than increased aerenchyma formation, we conclude that competition among roots is significant in the reference phenotype. Besides increasing root cortical aerenchyma formation, other anatomical phenotypes, such as reduced cortical cell file number and increased cortical cell size reduce root respiratory demands, potentially increasing carbon use efficiency even more if integrated with the SCD phenotype (Lynch, 2015).
Besides improving plant development under drought to mitigate the impact of climate change on yield, it has been suggested that improved carbon storage in the soil by plant roots would sequester atmospheric CO 2 (Lorenz and Lal, 2005;Dondini et al., 2009;Kell, 2011;Kell, 2012). Due to the reduced presence of oxygen and biological activity, carbon recalcitrance increases with soil depth (Lorenz and Lal, 2005). While OpenSimRoot does not yet simulate the carbon cycle in the soil, we can calculate the carbon deposited below 50 cm in the soil in the form of roots and exudates. In all environments, the SCD phenotype notably increased carbon deposition in deep soil. Increasing the carbon deposition by at least 0.3 gram per plant, as we observed in Zaria and Jalisco, at the planting density of 8 plants per m 2 amounts to deposition of 24 kg of carbon, or removal of approximately 88 kg of CO 2 per hectare. This increase of 0.3-gram carbon per plant deposited does not take into account respiration and is observed after 42 days of growth, so the actual amount after harvest (80 to 160 days) will be greater. Our simulations suggest that the SCD ideotype could increase shoot dry weight as well as carbon deposition in deep soil. It should be noted that increased carbon deposition in the soil does not necessarily mean more carbon will be sequestered long term, and carbon deposition can even lead to overall carbon losses from the soil (Fontaine et al., 2007). More detailed studies combining detailed plant root models and soil carbon cycles are needed to determine in what circumstances the SCD ideotype is beneficial for carbon sequestration.
Our results highlight the importance of soil-plant interactions when modeling crop development since we see major differences among the three soils. Soil compaction can negatively impact root growth and plant development (Taylor and Brar, 1991;Unger and Kaspar, 1994). In the Iowa soil, the soil bulk density below 45 cm is greater than in the other soils, which leads to large differences in root length in deep soil strata, as compared to the other two soils. Where the reference phenotype has at least 1000 cm and up to 3000 cm of root length below 50 cm in Zaria and Jalisco, in Iowa, it was less than 100 cm. Soil compaction caused by long term use of heavy agricultural machinery, as occurs in the Iowa environment, makes it more difficult for axial roots to reach deeper soil layers, and leads to a reduction of lateral root growth in the compaction layer, resulting in a more shallow root system. This is consistent with the results of a field study where 94% of total maize root length was located in the upper 60 cm of soil, with 84% in the top 30 cm, and water uptake mainly occurred in the upper 75 cm of soil (Laboski et al., 1998). Rather than the maximum rooting depth, the depth of the most densely rooted soil layer was important for maize plants to cope with drought (Yu et al., 2007). Irrespective of the presence of a compaction layer, the SCD phenotype increases subsoil root length notably (rather than just the maximum rooting depth), explaining the increased drought resistance.
The soil affects the way roots grow and in turn, roots affect soil properties. Plants growing in compacted soil will have a more shallow root system and the water they take up from the topsoil increases soil penetration resistance so that reaching deeper soil layers where water might still be available becomes more difficult (Colombi et al., 2018). Our results suggest that the SCD ideotype leads to greater water content near the roots than the reference phenotype at most depths, resulting in reduced soil penetration resistance. In drought environments, the reference phenotype has greater water content around the deepest roots than the SCD phenotype had around the same depth, but the SCD phenotype had similar water content in deeper soil strata (the SCD phenotype had deeper roots). This makes sense because root length density will be very low at the deepest point of a root system, so the plant will not take up much water there. Because the SCD phenotype is deeper than the reference phenotype, it has access to deeper soil water and water capture has greater vertical distribution. In this way, the SCD ideotype can help plants avoid the feedback loop (shallow root systems take up more water from the topsoil, which increases topsoil penetration resistance, making it harder for roots to reach the subsoil) described in (Colombi et al., 2018). Steeper axial roots mean that a root system is less wide but at commercial planting densities a wide root system will lead to increased competition with neighboring plants and is unlikely to increase water capture at the stand scale.
If plants with the SCD root phenotype perform better under drought without any tradeoffs one would expect maize plants to already possess this phenotype. That substantial variation exists for root phenotypes in maize suggests that a single root phenotype may not be optimal in all conditions. If water and nitrate are readily available then the benefit of the SCD ideotype for water and N acquisition is irrelevant, but the reduced carbon demands remain beneficial. Previous studies have shown that in soils with low phosphorus availability, high root length density in the topsoil is beneficial (Lynch and Brown, 2001;Rubio et al., 2003;Postma et al., 2014). While nitrate is mobile in the soil, if precipitation rates are low then nitrate will remain concentrated in the topsoil and shallow root systems perform better (Dathe et al., 2016). It has been suggested that plants with an SCD phenotype may also be more susceptible to root loss caused by diseases and herbivory, because each root represents a larger fraction of the total root system (Schäfer et al., 2022a).
Phenotypic root plasticity is the ability of a plant to alter its root phenotype in response to a change in environmental cues (Schneider et al., 2020). Broadly, phenotypic plasticity can be either apparent (i.e., those lacking adaptive value -e.g., reduced root biomass in response to stress) or adaptive (i.e., those associated with the plant fitness -e.g., increased lateral root density in response to local nitrogen availability) or even maladaptive (such as reduced root hair length under nitrogen stress). OpenSimRoot has the capability to capture both apparent and adaptive plasticity (as well as maladaptive plasticity if that is of interest). However, in this work, apparent plasticity (i.e., reduction in shoot and root growth reduces in response to drought) is captured while the adaptive plasticity is set to null (i.e., the predefined state of the root phenes does not actively change in response to droughte.g., nodal root number or lateral densities do not alter in response to drought). Note that because OpenSimRoot simulates resource dynamics, a phene state that alters the resource capture (in this case, water) will have feedbacks on root growth. Adaptive root plasticity is a spatiotemporally dynamic phenomenon (Schneider et al., 2020). In other words, adaptive root plasticity varies with root age, root type, and degree of the perceived environmental cues, which are spatially and temporally dynamic and much more complex. Certainly, root architecture shows adaptive plasticity in response to environmental changes, however, not all root phenes are plastic in nature. Being under genetic control, the plastic behavior of a root phene is often genotype-specific (i.e., a phene could be plastic to an environmental cue in one genotype but non-plastic in another genotype of the same species) (Schneider et al., 2020). Evaluation of such adaptive plasticity for each of the considered root phenes in this work would thus require a multitude of simulations. Such in silico exploration of root plasticity is possible but falls beyond the scope of this work and deserves to be a separate study.
The results presented here should be understood to be obtained with a number of limitations. As in any modelling study, many processes have been simplified or omitted. While care was taken to include factors relevant to modeling photosynthesis under drought such as canopy temperature, canopy self-shading, light scattering, and stomatal conductance in response to plant water status, the model is inevitably a simplification of reality. We assumed the stomata can close and open instantly and the canopy is always at the equilibrium temperature. OpenSimRoot does not simulate phytohormones or explicit nutrient and assimilate transport through the plant. While the effect of impedance on root growth is now taken into account, root growth is not affected by other environmental factors such as soil O 2 concentration or soil temperature, though we do not expect these factors to have a major impact on the results of this study. Due to computational constraints, we are simulating at finite temporal and spatial resolutions: the simulations run at a timestep of 0.1 days, which means we simulate the plant at different times of day (morning, noon, afternoon, etc.) while keeping computational costs manageable. The soil has a resolution of 1 cm, which allows for some differentiation between the resource capture of different roots but means the soil water gradient around individual roots cannot be accurately simulated. A limitation to the current model is that while the water flowing out of the roots during a hydraulic lift is simulated, this water is not added back to the soil, which means the topsoil has lower water content than it should. The total outflow is only a few percent compared to the total water content in the soil or the total change in soil water over the 42 simulated days so this will have a minor effect on results.
Perhaps the greatest limitation is that plants are only simulated for 42 days after emergence instead of until harvest. This is because the exponential increase in root segments over time means that computational costs of every simulated day increase as the root system grows. Simulating a maize plant for 42 days requires several days of computational time in some cases. While shoot dry weight correlates with yield (Hay, 1995;Huetsch and Schubert, 2017), the relationship between them is nonlinear, especially under drought stress, and there is no guarantee that the phenotypes with the greatest shoot dry weight at 42 days will also have the greatest shoot dry weight at harvest. Therefore, our results are relevant to vegetative growth under drought rather than yield under drought. An important factor under terminal drought is the banking of soil water to sustain reproductive growth (Zhang et al., 2019), a process not simulated here. Likewise, for carbon sequestration, the final amount of carbon deposited in the soil is relevant and it is possible a plant starts off slowly but deposits more carbon overall. However, since we see a strong correlation between shoot dry weight and carbon use efficiency under drought, the qualitative conclusions drawn from our results should be valid even when older plants are considered.
In summary, OpenSimRoot_v2 is capable of simulating maize growth under water deficit stress, encompassing root exploration of drying, hardening soil, water acquisition, and photosynthetic responses to water availability. Implementing this model, we present in silico evidence that more parsimonious root phenotypes, as epitomized by the SCD ideotype, that optimize the metabolic costs of soil exploration and water acquisition are capable of improving vegetative growth and carbon deposition in deep soil under drought in distinct maize production environments.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions
ES, IA, and JL contributed to conception and design of the study. ES and IA selected, implemented and tested the models, ran the simulations and analyzed the data. EF, LB, and MO supervised model implementation and helped resolve modelling issues. ES wrote the first draft of the manuscript. ES, IA, and JL wrote sections of the manuscript. All authors contributed to manuscript revision, read and approved the submitted version.