Simulating grass phenotypic plasticity as an emergent property of growth zone responses to carbon and nitrogen metabolites


 Phenotypic plasticity - the ability of one genotype to produce different phenotypes depending on growth conditions - is a core aspect of the interactions between plants and the environment. The model CN-Wheat simulates the functioning of a grass culm and the construction of traits as properties emerging from the feedback loops between morphogenesis, the environmental factors and source–sink activities. The plant is seen as a self-regulated system where leaf growth is driven by carbon and nitrogen metabolism within each leaf and by coordination rules between successive leaves. Here, we investigated the ability of this approach to simulate realistic grass phenotypic plasticity and explored plant behaviour in a wide range of growth conditions.The growth of grass monoculms, with traits similar to a wheat stem, was simulated for highly contrasting conditions of soil nitrogen concentration, incident light and planting density. The monoculms were kept vegetative and produced ~15 mature leaves at the end of the simulations. The model simulated highly contrasting phenotypes. Overall, the simulated trends and the magnitude of responses of leaf and plant traits to growth conditions were consistent with the literature on grass species. These results demonstrate that integrating plant functioning at organ scale can simulate, as an emergent property, the phenotypic plasticity of plants in contrasting light and nitrogen conditions. Besides, simulations of the internal variables of plants gave access to plant trophic status across plant ontogeny and plant environments. In conclusion, this framework is a significant step towards better integration of the genotype-environment interactions.


• Gauthier et al.
To date, most crop models are based on a common approach where carbon (C) acquisition is considered the limiting step in biomass accumulation (Poorter et al. 2013). A similar observation can be made for functional-structural plant models (FSPMs) as the potential growth rate of an organ (sink strength) is usually met only in case of sufficient resource availability (mostly C). In addition, key phenotypic traits, such as the phyllochron or specific leaf area, are usually simulated by empirical relations. The ability of models to account for the phenotypic plasticity is therefore highly related to the range of experimental conditions used to calibrate these relations. This approach is limiting our ability to explore the response of plants to the wide range of novel and upcoming growth conditions since they cannot be fully addressed experimentally . This calls for more mechanistic-and process-based models of the whole-plant functioning where plant and leaf traits are emergent properties resulting from the response of biological processes to the environmental factors. Functional-structural plant models simulate the behaviour of individual plants from an explicit description of their architecture which interacts with the environmental factors and plant functioning (Godin and Sinoquet 2005;Vos et al. 2010;Louarn and Song 2020). As such, FSPMs provide a suitable framework to explore the phenotypic responses of plants to a wide range of conditions.
In grasses, many FSPMs have been developed to simulate the growth and development of leaves in maize, rice, wheat or perennial grasses as the dynamics of leaf area expansion strongly determines plant behaviour and agronomic performances (Monteith 1977;Gallagher et al. 1979). The leaves of grasses have the particularity of growing from the base of the plant, sheltered from the light in a socalled pseudostem that is made of the sheaths of the older mature leaves (Kemp 1980b;Schnyder et al. 1987). Although the causal link remains to be established, several studies concluded that the emergence of leaf tips out of the pseudostem influences the ontogeny of the growth zone as well as variations in the leaf elongation rate of the emerging leaf or younger leaves. The first evidence supporting this hypothesis is that the pseudostem length strongly influences the final length of the leaves that growth within (Wilson and Laidlaw 1985;Casey et al. 1999). Secondly, many synchronisms have been found between the emergence of a leaf tip outside the pseudostem and the timing of transitions in the kinetics of elongation of younger leaves (Williams 1974;Malvoisin 1984;Tesařová et al. 1992;Skinner and Nelson 1995;Fournier 2000;Ljutovac 2002;Verdenal 2009). These observations led to the implementation of coordination rules in FSPMs in order to account for the effect of plant architecture on leaf morphogenesis, leading to the concept of self-regulated architecture (Durand et al. 1998;Fournier and Andrieu 1999;Fournier et al. 2005;Verdenal et al. 2008;Zhu et al. 2014;Gauthier et al. 2020;Vidal and Andrieu 2020). In addition to these synchronisms, C and N are the two main resources, together with water, that limit growth McIntyre 2001;Dornbusch et al. 2011). For grasses cultivated at contrasting N fertilization, the literature reports correlations between the concentrations in C and N in the leaf growth zone and the rate of leaf elongation (Kemp 1980a;Volenec and Nelson 1984;Gastal and Nelson 1994). Schnyder and Nelson (1989) also showed that tall fescue cultivated under low irradiance had smaller lamina width and specific leaf mass, which was associated with a lower carbohydrate concentration in the leaf growth zone. Therefore, it seems that processbased models would benefit from explicitly considering the amounts of C and N substrates available for leaf growth. This calls for models describing the C and N economy of the whole plant in interaction with the environmental conditions.
Recently, we developed CN-Wheat, a mechanistic model simulating the whole-plant functioning where plant and leaf traits are emergent properties originating from the feedbacks between the morphogenesis, local environmental factors and source-sink activities described at organ scale (Barillot et al. 2016a, b;Gauthier et al. 2020). In CN-Wheat, the plant is seen as a self-regulated system where leaf growth is driven by C and N concentration within each leaf and by coordination rules between successive leaves. Although CN-Wheat was calibrated for wheat, the representation of plant architecture, the existence of synchronisms between successive leaves as well as the modelling hypotheses about leaf growing zones and the modulation of their activity by substrate content are assumed to be generic to most grasses, particularly straw cereal crops and perennial grasses. In a previous study, we demonstrated the ability of CN-Wheat to simulate the behaviour of a wheat plant during the early vegetative stages in a usual agronomic condition (Gauthier et al. 2020). To further assess the model, the aims of the present work were (i) to evaluate the ability of CN-Wheat to simulate realistic dynamics of leaf and plant traits in response to contrasting growth conditions and (ii) to improve our knowledge on plant functioning under a wide range of growth conditions from the simulations of the internal variables of plants as provided by CN-Wheat. The present paper focuses on the plasticity of major foliar traits (e.g. dimensions, biomass, specific leaf area, N content). The mechanisms driving tiller emergence and regression are not yet implemented in CN-Wheat (although a preliminary and parametric version was used in Gauthier et al. 2020). Given the large diversity of growing conditions explored in the present study, we therefore chose to simulate leaf growth plasticity for monoculms of grass plants instead of introducing putative rules of tiller emergence. The consequences of this choice on the simulated plant plasticity are discussed later in the paper.

M ATER I A L S A ND M ETHODS
The growth of wheat plants was simulated with CN-Wheat in a wide range of conditions with differing planting densities, irradiances and soil N concentrations. Because CN-Wheat does not capture the mechanisms that define tiller emergence, we simulated wheat monoculms, thus focusing our analysis on the plasticity of leaf growth. The monoculms were kept vegetative, meaning that they produced an indeterminate number of phytomers and that internodes remained short. We present below: (i) the key aspects of CN-Wheat and the adaptations made to the initial version of the model as presented in Gauthier et al. (2020), (ii) the simulation conditions and (iii) the calculation of indices to compare simulations with observations. Gauthier et al. (2020) presented the dynamic version of CN-Wheat, which was calibrated to simulate the growth of a wheat plant under usual field conditions during the early vegetative stages. Preliminary simulations using the former version of CN-Wheat (Gauthier et al. 2020) revealed that some aspects of the model had to be adapted to Simulating grass phenotypic plasticity • 3 simulate the growth of monoculms that produce a high number of leaves under a wide range of environmental conditions as tested here. Below, we present an overview of CN-Wheat model and the processes that were evaluated in this study. Then, the changes to the model are briefly introduced, while a detailed description and underlying motivations are presented in Supporting Information-Data S1.

Model description
2.1.1 Overview of the model. CN-Wheat integrates C and N metabolism at organ level and simulates plant morphogenesis based on C and N metabolite concentrations. The main biological processes involved in the primary metabolism of C and N are described at organ level and are driven by the local concentrations of metabolites and the local environmental conditions (mainly light and temperature). Light distribution within the 3D canopy is calculated with a radiative model using the projective method (Chelle and Andrieu 1998). Briefly, CN-Wheat covers the processes of photosynthesis and transpiration, growth of each shoot organ in mass and dimensions and growth of the roots in mass, senescence and the fluxes of C and N within the plant (transport, synthesis/degradation, losses by respiration and root exudation). A summary of the model is given below regarding plant structure and shoot morphogenesis; the complete description of the model and its underlying hypotheses are available in Barillot et al. (2016a) for the C and N metabolism within a static architecture, and in Gauthier et al. (2020) for the extension of the model to a dynamic architecture.

Plant structure.
Each plant is represented as an ensemble of growing and mature shoot phytomers, a unique root compartment and a phloem compartment that is a common pool of mobile metabolites accessible to all organs of the plant. Each mature shoot phytomer has an internode, a sheath (with a distinction between the exposed and the enclosed tissues that face contrasting light environments) and a lamina. The growing leaf is represented by two types of compartments: the zone of the growing leaf that is hidden inside the pseudostem (hz) and the emerged tissues. During the simulation, new phytomers are initiated.
Each organ has several compartments describing their structural mass (mainly cell wall compounds) and the main metabolites involved in the primary metabolism of C and N. The model considers 'precursor' metabolites (triose phosphates, nitrates), mobile metabolites (sucrose, amino acids, hereafter denoted AA) and storage metabolites (fructans, starch, proteins). The shoot organs are also characterized by their dimensions and 3D shape, which are used to calculate light distribution, whereas the root compartment does not have an explicit architecture.

Leaf growth. The rate of phytomers initiation (plastochron)
is solely depending on the temperature of the shoot apical meristem. Leaf growth takes place in the hz, located inside the pseudostem, while the emerged tissues of the growing leaf are considered mature, photosynthetic and importers of root nitrates according to their transpiration rate. Sucrose and AA concentrations in the hz drive leaf elongation rate and modulate lamina width and the specific structural masses (SSMs, g m −2 ) of laminae and sheaths. In turn, leaf growth consumes sucrose and AA in the hz for structure synthesis and associated respiration.
Sucrose and AA are imported in the hz from the phloem and from the emerged tissues of the leaf (if any) following diffusion laws.
In addition, the model accounts for a self-regulation of the plant architecture, as the elongation of one leaf is coordinated with the elongation of the previous leaf. More specifically, the elongation of the whole leaf (lamina and sheath) is modelled as a two-phase kinetic with the transition from the first to the second phase being triggered by the emergence of the previous leaf tip outside the pseudostem. After initiation, leaf length follows an exponential-like function (phase I) which relative elongation rate (RER) is regulated by sucrose and AA concentrations in the hz. At the emergence of the previous leaf outside the pseudostem, leaf elongation enters a sigmoidal kinetic (phase II) which shape is mostly predefined. The duration of phase II only depends on temperature. The rate of leaf elongation during phase II depends on the length of the leaf at the beginning of phase II and is only slightly modulated by sucrose and AA concentrations in the hz, reflecting the fact that final leaf length is largely defined at previous leaf emergence. The final lengths of lamina and sheath are calculated according to an ontogenic sheath:lamina ratio based on leaf rank. We defined leaf ligation as the moment when the ligule emerges from the pseudostem, which corresponds to the complete emergence of the leaf blade.
From leaf initiation to leaf emergence, the leaf is only made of an hz whose structural mass is calculated from an empirical function of its length. At the emergence of the leaf, the SSM of the lamina and the sheath and the ratio between the lamina maximum width and its length (used for the leaf 3D reconstruction) are defined according to empirical functions of sucrose concentration in the hz averaged over the previous two phyllochrons. The structural masses of the emerged laminae and sheaths are then dynamically calculated from their respective SSMs and dimensions.

Modifications related to leaf and root growth.
First, the empirical functions that defined the lamina maximum width and structural specific lamina mass (SSLM) in the former version of the model were not intended to cover a large range of environmental conditions. We therefore modified the response laws of these leaf traits to sucrose concentration in order to have realistic results for extreme growth conditions [see Supporting Information- Table S1.1]. Second, we introduced an upper bound to leaf length prior to previous leaf emergence (phase I), which reflects a maximum length of the meristematic zone of the leaf [see Supporting Information- Table S1.2]. This change was necessary to simulate realistic leaf length in case of high number of leaves. Additional changes on leaf and root growth are detailed in Supporting Information-Data S1.

Modifications related to the photosynthesis.
Preliminary simulations using the initial calibration of the model (Gauthier et al. 2020) revealed that the photosynthesis rate was unlikely high for conditions of high irradiance and for conditions resulting in high specific leaf N (SLN, g N m −2 ). Consequently, the initial photosynthesis submodel was adapted by (i) implementing a negative feedback of non-structural carbohydrate (NSC) concentration on photosynthesis, (ii) changing the N dependences of the initial photosynthesis submodel. Such changes were possible because CN-Wheat simulates the concentration 4 • Gauthier et al.
of metabolites in each organ, which allows implementing detailed regulations that are not considered in most models. Specifications of the changes to the initial photosynthesis submodel are given in Supporting Information-Data S1 and briefly presented below.

Virtual experiments
The model parameters (~160) were kept identical to those presented in Gauthier et al. (2020), with the exception of the changes that have been described in the previous section. During the simulation, no floral transition was simulated meaning that the culm continuously produced new vegetative phytomers and that the internodes stayed short.
As in Gauthier et al. (2020), the model was run with an hourly time step, except for the distribution of light between phytoelements, which was updated every 4 h to save computational time. This means that the absolute absorption of photosynthetically active radiation (PAR) by each organ was calculated hourly, using hourly values of incident PAR, but that the relative distribution of PAR within the 3D canopy was considered constant during a 4-h period. Each simulation run in ~2.5 h on a personal computer. In total, the virtual experiments consisted in 80 simulations and ~500 preliminary simulations were performed. Thus, the model has been deployed on a High-Performance Computing centre (MESO@LR-Platform, Université de Montpellier, France), which allowed to perform all the simulations in parallel.
Below, we present the environmental conditions of the simulations and the initialization of the simulated plants. Indices characterizing plant growth and trophic status were calculated from simulation results and are present in Supporting Information-Data S1.

Environmental conditions CN-
Wheat was run to simulate the growth of a monoculm for 3500 h (~145 days), the day/night air temperatures were fixed at 13/7 °C and the photoperiod at 12 h. Soil temperature was constant at 10 °C. Atmospheric CO 2 concentration, relative humidity and wind speed at the top of the canopy were set at 360 ppm, 75 % and 3 m s −1 , respectively.
The main virtual experiment (VE1) consisted in 84 simulations with three planting densities, six soil N concentrations and seven irradiances (Table 1). Planting densities were 1, 500 and 750 culms per m 2 . For simplicity, the density of 1 culm per m 2 will be referred below as isolated plants. For comparison, a commercial wheat field in France would typically count between 500 and 600 culms per m 2 at flowering. The NO 3 − concentration of the soil was set at constant values (thus mimicking hydroponic conditions), ranging from 0.5 to 10 g m −3 . Since the model considers a soil of 1-m depth, such concentrations are equivalent to constant soil N concentrations ranging from 5 to 100 kg ha −1 . For comparison, 80 kg ha −1 is a concentration that can be found immediately following fertilizations in a commercial wheat field. However, in such field, the concentration would decrease rapidly with plant uptake (Devienne- Barret et al. 2000). Finally, four incident PAR levels were simulated, from 3 to 11 mol m −2 day −1 (80-250 µmol m −2 s −1 ) for isolated plants and from 8 to 65 mol m −2 day −1 for plants in stands (175-1500 µmol m −2 s −1 ). Different ranges of incident PAR were applied to isolated plants and plants in stands to account for the differences in light interception per plant. For comparison, PAR of 22 mol m −2 day −1 with 12-h photoperiod and 10 °C daily mean air temperature corresponds to the usual weather conditions in Grignon (France) in March.
In addition, we performed another virtual experiment (VE2) to assess the effect of the initial state of the model on the sensitivity of the simulations. VE2 consisted in eight simulations under a subset of environmental conditions taken from VE1 (two planting densities, two soil N concentrations and three irradiances; Table 1).

Model initialization
In VE1, the initial state of the main stem and roots (Init VE1 ) was defined as in Gauthier et al. (2020) in terms of dimensions, masses and metabolite concentration Briefly, this initial state reflects the state of wheat plants grown in the field in Grignon at the stage of leaf 4 emergence (mid-December). In VE2, a second initial state of the plant Init VE2 was defined by decreasing the initial shoot N fraction to 2.4 % (compared to 4.4 % in Init VE1 ) while keeping the same repartition of N within the plant as in Init VE1 .

R E SULTS
We focus on the simulation results of VE1, which corresponds to the growth of vegetative monoculms under 84 conditions contrasted by planting density, soil N concentration and incident PAR. The simulation results of VE2 are partially presented to illustrate the effect of the initial plant state on internal concentrations. For clarity's sake, we refer below to growth condition treatments by short names, e.g. D1_N5_ PAR4 refers to a density of 1 plant per m 2 , a soil N concentration of 5 g m −3 and an incident PAR of 4 mol m −2 day −1 .
First, we present the results of plant production in terms of biomass and leaf area. Then, we give an outlook on C and N concentrations in the phloem, which reflect C and N availability at plant scale. Finally, we present the phenotypic plasticity at leaf scale and then at plant scale.

Plant production
At the end of the simulation, plants reached highly contrasting sizes and masses depending on growth conditions (Fig. 1). Total area of laminae per plant ranged from 20 to 325 cm 2 at the end of the simulation, while the living dry mass ranged from 0.2 to 7.1 g. Both total Simulating grass phenotypic plasticity • 5 lamina area and plant dry mass were higher with increasing incident PAR, and with increasing soil N concentration. The effect of soil N concentration was mainly observed up to 4.5 g m −3 and then saturated for higher concentrations (Fig. 1). Isolated plants and plants in stands showed some differences in their response to the growth conditions. A 50 % change in total lamina area was associated to a 50 % change in dry mass for plants in stands, but to a higher change, up to 70 % for isolated plants (Fig. 2). Thus, plant dry mass was more sensitive to the variations of lamina area induced by environmental conditions in isolated plants than in stands (Fig. 2), reflecting that in stands, differences in lamina area did not imply additional differences in intercepted PAR because of canopy closure.
However, when considering the phenotype of the plants at the end of the simulations, the effect of increasing plant density closely resembles that of decreasing incident PAR, which is consistent with experimental observations (Poorter et al. 2019). Thus, in the following sections, we present the simulation results for D500 only. In addition, variations with soil N concentration were monotonous. So, for clarity' sake, when presenting variables across time or per phytomer rank, the effect of soil N concentration is exemplified for N0.5 and N5, that are very low and very high soil N concentrations, respectively.

Internal concentrations
In the context of simulating phenotypic plasticity, the concentrations in C and N of the hz are of special interest since, in our model, sucrose and AA concentrations drive the leaf elongation rate, lamina width and specific mass of lamina and sheath. Also, the concentrations of sucrose and AA in the phloem are an indicator of the trophic status of the plant as they represent the amounts of mobile C and N that can flow from sources to sinks. In our simulations, the concentrations of sucrose and AA in the hz were always in equilibrium with the phloem because of low flux resistance between these compartments. Below, we exemplify the dynamics of sucrose and AA concentrations of the phloem during the simulation for D500 at four PAR × N conditions.

Sucrose. 3.2.1.1 Ontogenic pattern.
Sucrose concentration increased at the beginning of the simulations, with a marked peak at ~500 h (leaf 6 emergence), ranging from 400 to 3400 µmol g −1 depending on the treatment ( Fig. 3A and B). Then, the growth of leaves 6 and 7 consumed sucrose in the hz (data not shown), which caused a decrease of the sucrose concentration until 1500 h in all treatments. From 1500 h onward, sucrose concentration of plants at D500 tended to stabilize between 200 and 800 µmol g −1 depending on the treatment ( Fig. 3A and B). In contrast, sucrose concentration of isolated plants continued to increase and reached similar values at 3500 h to those at 500 h (not shown).

Effect of PAR.
The sucrose concentration increased with incident PAR ( Fig. 3A and B). Concentration differences between PAR treatments started from the beginning and were maintained until the end of the simulation. At low PAR, fluctuations in sucrose concentration resulted in transient very low values at ~1200 h and ~2500 h.

Effect of N.
Overall, the effect of soil N concentration on sucrose concentration in the phloem was weaker than that of incident PAR. Figure 3A and B shows that sucrose concentration started to be higher with increasing soil N concentration from ~1000 h.

Effect of initial state of the plant.
The initial state of the plant strongly impacted the sucrose concentration of the phloem during at least the first 800 h of the simulations, with a marked effect on the initial peak. Differences no longer existed in high PAR treatments after 800 h, while differences existed until 2000 h at low PAR and 2500 h in conditions of low PAR-low N.

Ontogenic pattern.
As for sucrose concentration, AA concentration showed a peak at ~500 h, especially for high incident PAR and high soil N concentrations ( Fig. 3C and D). After 1500 h, AA concentration increased in situations of high incident PAR and high soil N concentration, e.g. reaching 600 µmol N g −1 at 3500 h in the simulation D500_N5_PAR1000 (Fig. 3D). The strong accumulations of AA after t = 1500 h were found to be associated with a rapid growth of the roots, which was due to high sucrose concentrations in the phloem and high allocation to the roots. In the other cases, AA concentration was relatively stable from 1500 h onward. The local minimums in AA concentration (and also in the total amount of AA per plant-not shown) appeared to be synchronized with the ligulation of successive leaves ( Fig. 3C and D). It seems that the transient increases were due to the ligulation of young leaves, which started to load AA towards the phloem. Then, the following decrease in AA concentration was probably due to the rapid growth leaf n + 1 (soon after leaf n ligulation), which consumed AA for the synthesis of structural mass and of proteins.

Effects of PAR and N.
Amino acid concentration increased with incident PAR and soil N concentrations ( Fig. 3C and D) from 150 h. The plant N fraction also increased with increasing incident PAR [see Supporting Information- Fig. S2.1], indicating that additional PAR led to a higher increase in N absorption than in C assimilation.  Simulating grass phenotypic plasticity • 7 Such simulation result departs from experimental results of Gastal and Saugier (1986) who observed a decrease in the shoot N fraction with increasing incident PAR for tall fescue cultivated in hydroponics. Although their N nutrition regimes are comparable to those applied in the present study, the former results were obtained in a different species and not for monoculms, which probably restricted the aboveground development in or simulations. These seemingly opposing results may therefore instead highlight that the effect of PAR on plant N content varies considerably depending on the respective responses of above-ground (especially tillering) and below-ground development in response to C (and possibly N) concentrations.

Leaf length. 3.3.1.1 Ontogenic pattern.
Overall, the pattern of leaf length along the stem was the same for all growth conditions, final leaf length generally increased with phytomer rank up to leaf 12 and was then almost stable (Fig. 4). Exceptions were observed for some phytomer ranks, e.g. leaf 10 was always shorter than leaf 9. Such exceptions contrast with experimental observations as smooth, monotonous leaf profiles are observed in the absence of transient stress for the leaves beard on internodes that remain short (Williams 1960;Gallagher 1979;Dornbusch et al. 2011;Abichou et al. 2018). The relative stability of leaf length for phytomers 12 and above, which was simulated in most growth conditions revealed that, before previous leaf emergence, the hz reached the upper bound allowed in the model. Such stability of leaf length for high phytomer ranks is often observed on perennial grasses (Verdenal 2009). Length of leaves 1-6 was almost insensitive to growth conditions as their potential final lengths were already set at model initialization.

Effect of PAR.
Leaf length increased with incident PAR and reached a maximum of 54 and 60 cm for the highest leaf ranks in D500_N0.5 and D500_N5, respectively. For high PAR values (above a threshold that depends on density and soil N concentration), leaf length was almost stable with the incident PAR. The simulated increase of leaf length with increasing incident PAR resulted from the upregulation of leaf elongation rate by sucrose concentration mainly, and to a lesser extent by AA concentration. This effect however saturated above a PAR threshold because sucrose and AA concentrations in the hz reached levels saturating the response of the leaf elongation rate. Literature reports a small positive effect of PAR on leaf length (Friend et al. 1962;Ljutovac 2002;Dornbusch et al. 2011), except for very low light levels which increase leaf length, e.g. by Friend et al. (1962) on wheat, by Allard et al. (1991) on tall fescue, by Paciullo et al. (2017) on Guinea grass and by Bos et al. (2000) on maize. This aspect will be addressed in the discussion.

Effect of N.
High soil N concentration strongly increased leaf length as extensively reported in the literature (Robson and Parsons 1978;Wilman and Mohamed 1980;Volenec and Nelson 1983;Fricke et al. 1997). Since we simulated constant soil N concentrations, it is difficult to compare quantitatively the magnitude of the simulated effects with the literature. However, on vegetative tall fescue and ryegrass, Wilman and Mohamed (1980) reported a 25 % increase of leaf length with N fertilization, which is of same order of magnitude than the length increase simulated here between N0.5 and N5 for intermediate leaves (rank 6-12) at D500_PAR1500.

Lamina maximum width and lamina mass per area.
Simulated lamina maximum width varied in the range of 0.3-3.0 cm, and lamina mass per area (LMA) varied in the range of 10-95 g m −2 . The highest simulated values of lamina maximum width and LMA are higher than those reported in experimental studies (Rawson et al. 1987;Dornbusch et al. 2011;Abichou 2016). In particular, we did not find experimental lamina maximum width values above 2 cm for wheat, nor experimental LMA values above 60 g m −2 for wheat and ryegrass. However, the simulations in the present work covered a larger range of growth conditions than those reported in the literature.

Ontogenic pattern.
For some growth conditions, the simulated values of lamina width and LMA presented sharp changes with phytomer rank [see Supporting Information- Fig. S2.2] due to the large ontogenic fluctuations of sucrose concentration ( Fig. 3A and B). This could mean that the model does not sufficiently buffer sucrose fluctuations and/or that the method used to calculate maximum lamina width and LMA is excessively sensitive to sucrose fluctuations.

Effects of PAR and N.
For all phytomer ranks, lamina maximum width and LMA markedly increased with incident PAR, while the effect of soil N concentration was lower; we illustrated the effects of incident PAR and soil N concentration on lamina 8, which presented intermediate responses ( Fig. 5A and C). As shown before, the incident PAR and soil N concentration affected sucrose concentration in the hz, which regulates the lamina maximum width and LMA. The magnitude of responses of the lamina maximum width and LMA varied markedly with the phytomer rank since sucrose concentration in the hz sometimes reached saturating levels of the response laws of ratio of lamina maximum width to potential final length [see Supporting Information- Fig. S1.1A] and of SSLM [see Supporting Information- Fig. S1.1B]. Thus, according to the model, organs that grow during a high source:sink (high sucrose concentration) are almost insensitive to growth conditions. The effects of incident PAR are consistent with the literature (Friend et al. 1962;Rawson et al. 1987;Bos et al. 2000;Dornbusch et al. 2011;Poorter et al. 2019). Friend et al. (1962) cultivated stands of spring wheat at the same mean daily temperature as our simulations, under a range of incident PAR. They observed that incident PAR had almost a linear and positive effect on the lamina maximum width (Fig. 5B) and LMA (Fig. 5D), which diminished for PAR above 22 mol m −2 day −1 . Overall, the simulations obtained from CN-Wheat were very close to these observations ( Fig. 5B and D). The slight positive response of lamina maximum width to soil N concentration (Fig. 5A) is consistent with the literature (Wilman and Mohamed 1980). In contrast, LMA tended to increase with increasing soil N concentration in our simulations (Fig. 5C) while the literature reports contrasting effects of soil N concentration on LMA depending on the grass species: a positive effect was observed in barley (Fricke et al. 1997) and wheat (Dreccer et al. 2000) but the opposite response was observed for rye grass (Robson and Parsons 1978) and tall fescue (Wilman and Mohamed 1980). Variables are expressed relatively to their value at PAR 22 mol m −2 day −1 . Dots are derived from experimental data for the third leaf of spring wheat seedlings cultivated on stands in growth chambers at 10 °C with a 24-h photoperiod (Friend et al. 1962).

Ontogenic pattern.
For all treatments, leaves emerged at a constant rate until leaf 9-10 and then at a slightly higher rate for the upper leaves, as exemplified in Fig. 6A. The rate of leaf ligulation followed the rate of leaf emergence. The stability of the time interval between the emergence (or ligulation) of successive leaves is an emergent property of the model, which appeared in all growth conditions. When averaged over the whole simulation period, the phyllochron ranged from 106 to 115 °Cd depending on growth conditions (Fig. 6B).

Effect of PAR.
The phyllochron decreased with increasing incident PAR and then stabilized above 43 mol m −2 day −1 of PAR. The negative effect of PAR on phyllochron is consistent with the literature, although the magnitude of the simulated response is smaller than that reported in the literature. Indeed, Triboï and Ntonga (1993) observed a 25 % increase in wheat phyllochron grown at low density in a greenhouse with 60 % shading. Likewise, Baumont et al. (2019) observed a 20 % increase in the phyllochron after a 40 % reduction of the incident PAR (16 vs. 10 mol m −2 day −1 ) for wheat plants cultivated at 28/24 °C with 16-h photoperiod. In comparison, the maximum response of the phyllochron simulated with the model consisted of a 10 % increase (at D500) after an 80 % reduction in incident PAR (65 vs. 11 mol m −2 day −1 ).

Effect of N.
The phyllochron also decreased sharply with increasing soil N concentration up to ~2.5 g m −3 and then stabilized. The literature reports an effect of soil N concentration on the phyllochron for high N deprivation only (Longnecker and Robson 1994). For instance, Longnecker et al. (1993) applied weekly doses of nutritive solution of several soil N concentrations on two wheat genotypes. They reported that the phyllochron increased by 5 % when reducing the soil N concentration from 22.4 to 4.2 g m −3 , and by ~15 % when reducing the soil N concentration from 11.2 to 0.7 g m −3 . The maximum simulated increase in the phyllochron was 4 % with a reduction of the soil N concentration from 10 to 0.5 g m −3 at D500. The latter concentration of 0.5 g m −3 led to strong N deprivation (low SLN, low shoot N fraction). We conclude that, even if a direct comparison with the literature is not possible, the simulated effect of light and N on the phyllochron might be too small although the simulated trends were correct.

Shoot:root dry mass ratio.
At the end of the simulation, the shoot:root ratio ranged from 2 to 6.2 at D500 (Fig. 7). Experimental shoot:root ratio from the literature on vegetative ryegrass during the first year of growth ranged from 1.5 to 8 for several treatments of N and PAR (Hunt 1975;Robson and Parsons 1978). Experiments on rice and wheat showed that manual detillering or a genetic inhibition of tillering resulted in a lower shoot:root (Hendriks et al. 2016;Nada and Abogadallah 2016). Such results must be considered when comparing simulated shoot:root values with the literature. Overall, the simulated values of shoot:root ratio seem consistent with the observations.

Effects of PAR and N.
Overall, the final shoot:root ratio decreased with incident PAR (above 22 mol m −2 day −1 ) and increased with soil N concentrations. These trends are in line with the literature (Hunt 1975;Robson and Parsons 1978;Poorter et al. 2012) and with the optimal partitioning theory, which states that plants allocate preferentially resources to the most limiting component, which is the shoot in case of C limitation and the root in case of N limitation (Bloom et al. 1985). In our simulations, this behaviour resulted from a limitation of shoot growth instead of being solely the response of shoot and root growth to C and N concentrations. Indeed, at the end of the simulation, leaf length reached maximum values in most growth conditions (Fig. 4), whereas root growth had no limitation other than sucrose concentration in the roots. In contrast, the final shoot:root ratio increased with incident PAR from 8 to 22 mol m −2 day −1 (Fig. 7), which corresponded to simulated plants for which leaf length did not reached maximum values (Fig. 4).

Radiation use efficiency.
Simulated radiation use efficiency (RUE), calculated from absorbed PAR and plant dry mass (including structural mass of the senesced tissues), varied between 1.5 and 3.5 g MJ −1 for D500. These simulated values are in the range of observed RUE and that used in some crop models (Sinclair and Amir 1992). In particular, the highest simulated RUE is consistent with the maximum RUE reported in the literature: 3 g MJ −1 for wheat and barley before heading (Gallagher and Biscoe 1978;Reynolds et al. 2012) and 4.2 g MJ −1 for vegetative ryegrass (Akmal and Janssens 2004).

Effects of PAR and N.
Simulated RUE increased with decreasing incident PAR and with increasing soil N concentration (mainly up to ~2.5 g m −3 ; Fig. 8), both trends are coherent with experimental observations (Sinclair and Horie 1989;Bélanger et al. 1992;Gómez et al. 2013). The model simulated a decrease of RUE with increasing PAR as a result of (i) the non-linear response of photosynthesis and (ii) of the accumulation of NSC with PAR, which further inhibit photosynthesis.

DISCUSSION
CN-Wheat relies on a bottom-up approach in which the growth processes are driven at a local scale by (i) the C and N metabolite concentrations of each organ, and (ii) coordination rules that set the timing of elongation of successive leaves as well as the timing of definition of key leaf traits (lamina width and SSMs of lamina and sheath). The integration of these local processes builds whole-plant traits (e.g. shoot:root ratio), temporal patterns (e.g. phyllochron stability) and spatio-temporal patterns (e.g. patterns of leaf dimensions along the axis). Here we explored the model's behaviour when simulating growth dynamics in contrasting environments and tested its ability to simulate the phenotypic plasticity of traits at leaf and plant scale.
To this purpose, the growth of vegetative monoculms cultivated under hydroponic-like conditions was simulated in a wide range of planting densities, soil N concentrations and incident PAR. The chosen growth conditions facilitated the interpretation of simulation results (N availability does not depend on planting density), and presented extremely contrasting growth conditions to experience the model. A drawback is that they did not allow simple comparison with literature that usually reports field experiments with plants producing tillers.

The model simulates realistic phenotypic plasticity for contrasting light availability and N fertilization
Overall, CN-Wheat simulated realistic phenotypic plasticity of leaf and plant traits. As a response to increasing incident PAR, our model simulated, at the whole-plant level, a decreasing fraction of dry mass allocated to roots and faster leaf emergence, and, at leaf level, thicker and wider leaves. The response to increasing planting density was largely similar to that of decreasing incident PAR, which was expected as soil N concentration was forced independently of plant development in our simulations. As a response to increasing soil N concentration, our model simulated, at the whole-plant level, a decreasing shoot:root dry mass and faster leaf emergence, and, at leaf level, increasing thickness, width and length. CN-Wheat simulates most leaf and plant traits as emerging properties of feedback loops between C-N acquisition and responses to metabolite concentrations, whereas most models estimate these key traits (e.g. RUE, shoot N content, LMA) from empirical relations. In addition, the present work demonstrates that the integration of these feedback loops at organ scale builds realistic responses to growth conditions, establishing an alternative to current models. We exemplify this property with two key traits at plant level that are usually set as inputs in most models. One example is the phyllochron: its stability across plant ontogeny in each growth condition emerged from the coordination rules that set the timing of elongation of the successive leaves, while simulated variations of the phyllochron to light  Simulating grass phenotypic plasticity • 11 availability and soil N concentration were coherent with the literature (Fig. 6B). Another example is the realistic RUE values and RUE variations (Fig. 8) that resulted from the balance between C and N assimilation and losses (respiration, tissue death and root exudation). Other examples include the shoot:root ratio, the pattern of leaf length along the axis, etc.
Three discrepancies between the simulated responses and the literature were identified. First, in response to low PAR, the model simulated shorter leaves (Fig. 4), whereas the literature frequently reports an increase in leaf length mainly due to the concomitant decrease in blue light irradiance (Gautier and Varlet-Grancher 1996;Christophe et al. 2006;Barillot et al. 2021). This highlights that this phenotypic response cannot be simulated from a model of plant morphogenesis regulated by C and N availability only, but photomorphogenetic processes must be implemented in the model. Second, phyllochron variations with environmental conditions were weaker than expected (Fig. 6B). Our coordination rules state that, at the emergence of the previous leaf, the leaf starts a sigmoidal elongation phase of constant elongation duration, whereas some experimental works reported variations of the elongation duration with light availability and N, suggesting that our model of leaf elongation might be further refined (Dale 1982;Fricke et al. 1997;Fournier and Andrieu 1999;Louarn et al. 2010). The absence of photomorphogenetic processes in the model probably also explains the limited response of the phyllochron to the low-PAR conditions (Gautier and Varlet-Grancher 1996). Indeed, sheaths generally become longer under low-PAR conditions, increasing the length of the pseudostem and thus the phyllochron. Third, the model simulated increasing plant N fraction with increasing PAR [see Supporting Information- Fig. S2.1], whereas the literature on hydroponic conditions showed opposite trend (Gastal and Saugier 1986). This result points out the need to revise some of the model formalisms related to N uptake and assimilation. It is nevertheless worth noting that, discrepancies between simulations and observations also likely originated from the very specific conditions of our simulations (mainly monoculm and hydroponic conditions).
Indeed, cultivating standard wheat cultivars under such contrasting conditions would have led to highly contrasting tillering levels (Friend 1965;Casal et al. 1986;Longnecker et al. 1993;Nelson 2000;Evers et al. 2006). Accounting for tillering in the model would likely lead to at least two main changes. First, tillers would have progressively increased the density of culms and the competition for light, leading to a different dynamic of the canopy structure compared to that produced by increasing the density of monoculms. Second, tillering plants are expected to have higher shoot:root ratio than plants with few or no tillers (Hendriks et al. 2016;Nada and Abogadallah 2016), which would likely result in a lower N:C ratio of the plant. Implementing the underlying mechanisms of tillering would allow further evaluation of the model against experimental data.

Changes in leaf traits and cyclic variations of mobile metabolite across plant ontogeny
Changes in the properties between successive phytomers reflect a combination from changes in the environment perceived by plant organs (e.g. light quantity and quality with closure of the canopy), changes in the internal state of the plant during the period when the properties of the phytomer are build (e.g. concentration of C and N, size of the sheath tube) and possible intrinsic changes in phytomer properties (heteroblasty). In the case considered here, a culm remaining vegetative, it is tempting to make the hypothesis that-apart from embryonic leaves for which additional complexity may be expected-differences between successive phytomers reflect a plastic response to a changing external and internal environment. In our simulations, most of the changes between successive phytomers are emerging properties of the response to evolving internal and external environment; the ability of the model to reproduce ontogenetic trends in traits is thus narrowly related to its ability to account for plant ontogenetic plasticity.
In most growth conditions, leaf length increased with phytomer rank, up to a plateau that depended on growth conditions. Noteworthy, the plateau did not show in the initial version of the model but was produced by introducing an upper bound to the size of the growth zone. This model parameter was necessary to reflect the fact that leaf lengths rarely exceed ~45-50 cm (with ~12 cm sheath) in several grass species like. wheat, barley, rye grass. The requirement for introducing explicitly a limit to the size of the growth zone likely results from constraints not taken into account in the model, such as the resistance to transport of water and nutrients through the immature growing zone (Verdenal 2009). For vegetative grasses that produce a large number of phytomers, experimental evidence is that the rank of the longer leaf can be found in various positions depending of growth conditions (e.g. 6 in Ljutovac (2002) and 12 in Robin (2011)). A decrease of length is frequently reported for ranks above this position, which may result from a decrease of light and N resources with closure of the canopy although alternative hypotheses may be speculated. Robin (2011) reported experiments in which vegetative ryegrass plants were grown up to 15-17 leaves in hydroponic conditions, with very low planting density and manual detillering that minimized shading. In these conditions, he observed a plateau of leaf length and width reached approximately at rank 12. This supports that, when plant development does not result in increasing competition for resources, leaf length increases with successive phytomers until reaching a plateau. Similarly in such conditions, both simulated and measured leaf width increase up to reaching a plateau.
Our model simulated a range of patterns for the change in metabolite concentration along time, with long-term impact of the initial plant state. For sucrose concentration ( Fig. 3A and B), long-term trend (at 500 plants per m 2 ) was a decrease across development, quite fast during the period of canopy closure, then slightly decreasing, as canopy intercepted all light, and canopy structural mass increased slowly as a result of growth and senescence. Superimposed with the long-term trend, a most unexpected pattern was some kind of cycling of sucrose concentration of decreasing amplitude. Initial N fraction of the canopy (defined from field experiments) was distant from the one that the simulated growth conditions would have produced. In response to that, the model did not describe a smooth trend towards equilibrium but rather a cyclic behaviour of decreasing amplitude and long, irregular period (several phyllochrons). Thus, even in stable conditions, the model showed a marked oscillating behaviour that resulted from the feedback loop between concentration and the dynamics of response to this concentration. This behaviour is what creates, in the model, differences in the magnitude of response of successive phytomers to environmental conditions and an irregular distribution of leaf width and mass under condition of low resources [see Supporting Information- Fig.  S1.1]. We have no precise data allowing to evaluate the fluctuations in sucrose but reported change in leaf traits in response to stress are generally relatively smooth, suggesting that the model overestimated these fluctuations. Indeed, the response functions of the construction of leaf traits in response to sucrose concentration were defined in a very pragmatic way, with probably too simplistic choice in the way to integrate over time the responses of leaf width and thickness to sucrose concentration. We preferred to keep these simple functions, which only intend to be preliminary proposals, because we lack precise information about (i) the mechanisms that drive the construction of leaf width and thickness, (ii) empirical data that link leaf width and thickness with internal plant parameters instead of environmental conditions. We expect that a mechanistic implementation of the processes that build leaf width and thickness would not lead to such oscillating behaviour as it would integrate the trophic status of the plant for a longer period.
For AA concentration, long-term trends depended on simulation conditions. In most conditions, AA concentration showed little longterm trend but combined high PAR and high soil N concentration resulted in significant accumulation of AA over time. The values simulated and the contrast depending on growth condition is consistent with little existing bibliographic data (e.g. Lawlor et al. 1987). Interestingly, AA showed marked cyclic oscillations, which were precisely coordinated with leaf ligulation. In most growth conditions, the amplitude of the oscillations represented a large share of AA concentration, meaning that the plant oscillated at one-phyllochron intervals between phases of high and low AA availability. Oscillations of N fraction of one-phyllochron intervals has been reported by Robin (2011) in barley plants grown in hydroponic. Robin (2011) investigated total N; thus, the relative variations were much lower, but assuming they came from AA, which represent a small fraction of total N, their amplitude was similar to that simulated here. Moreover, Robin (2011) mentionned that Irving and Matthew (unpublished results) also observed such oscillations in AA concentrations which may represent a timing mechanism for leaf developmental processes. In CN-Wheat, the timing of leaf development is driven by a speculative coordination rule derived from experimental evidence showing synchronisation between the leaf elongation phase and leaf emergence. It has been proposed that this synchronisation between leaves comes from a light or gas signal perceived at the time of emergence of a leaf (Casey et al. 1999). Our simulation and above-mentioned experimental results suggest that oscillation in AA concentration is an alternative or re-enforcing mechanism contributing to produce a stable rate of leaf appearance.

Plant's C and N economy is an emergent property of CN-Wheat and drives plant morphogenesis
Plant's C and N economy is a key aspect of CN-Wheat as it permanently interacts with plant morphogenesis. Many models simulate morphogenesis from a source:sink ratio, which is not measurable and considers only C (Luquet et al. 2006;Letort et al. 2008;Mathieu et al. 2009;Evers et al. 2010). Here, morphogenesis is driven by C and N metabolite concentrations at organ scale, which can be compared to experimental data.
Such approach implies that the model must duly account for the main metabolic processes that control plant's C and N economy, which are the processes of C and N acquisition, as well as C and N losses (respiration, exudation, tissue death), and C and N partitioning among the different metabolic compartments (different forms of metabolites and organ structure) and organs.
This work allowed identifying that the classic biochemical model of photosynthesis (Farquhar et al. 1980;FCB model) led to probably unrealistically high fractions of NSC and RUE [see Supporting Information-Data S1] in some extreme conditions when the model simulated very thick leaves with a very high specific leaf N (SLN). So, we switched from dependency of the FCB parameters on SLN to dependency on non-structural proteins, and we implemented an empirical feedback regulation by the NSC concentration. These implementations remain putative but are original as compared to other models based on the FCB model, and were possible since CN-Wheat simulates the concentrations of C and N metabolites at organ scale. The implemented feedback regulation of photosynthesis relied on an empirical function that decreases the rate of photosynthesis in function of NSC concentration, which has been calibrated from experimental results on wheat (Azcón-Bieto 1983). Because the underlying mechanism is not yet elucidated (Sharkey 2019;Dingkuhn et al. 2020), a mechanistic implementation is not yet considered.
Other processes have a strong impact on plant's C and N economy but have been less studied (e.g. leaf and root senescence). In addition, processes must present genotypic variations that result in a wide range of plant behaviour. For exploratory purposes, we performed additional simulations with a lower maximum rate of root growth or with the implementation of a limited root lifespan (not shown), that are little known aspects of plant morphogenesis. These simulations resulted in higher shoot:root ratio, lower AA concentration and lower plant N fraction and thus, revealed that the model behaviour is largely determines by parameters that modulate plant architecture and organs lifespan. Such parameters might play a key role in the differences of behaviour that are observed between species and cultivars, following the assumption that model parameters do represent genotypic characteristics (e.g. Dingkuhn et al. 2005;Luquet et al. 2012;.

CON CLUSION
CN-Wheat provides an original and explicit description of the processes underlying shoot morphogenesis, which provides new possibilities to link the phenotypic plasticity of plants to their C and N metabolic status at the place and time where traits are built. This unique level of integration makes CN-Wheat a good tool to explore plant functioning in contrasting environments and opens up new possibilities to account for genotype-environment interactions.

SUPPORTIN G INFOR M ATION
The following additional information is available in the online version of this article-Data S1. Additional information regarding materials and methods. Data S2. Additional results.