Jean-Franc ¸ Ois Soussana
email: jfsoussana@clermont.inra.fr

Vincent Maire

Nicolas Gross

Bruno Bachelet

Loic Pagès

Raphaël Martin

David Hill

Christian Wirth

Gemini: A grassland model simulating the role of plant traits for community dynamics and ecosystem functioning. Parameterization and evaluation

Keywords:

A structure-function-diversity model of grassland ecosystems (Gemini) has been developed. For a potentially unlimited number of clonal plant populations, it explicitly simulates competition for two key resources (light and nitrogen) along vertical canopy and soil profiles. Population turnover, shoot and root morphogenesis, photosynthesis, respiration, transpiration, N acquisition by uptake, allocation of assimilates between structural compartments, and reserve storage and remobilization, are simulated for each plant population. The object-oriented structure of the modeling framework allows to couple, or not, the simulated plant populations to other sub-models describing climate variables, soil functioning, grazing behavior and grassland management. Partitioning of growth between shoot structures, leaf photosynthetic proteins and roots is based on two assumptions: (i) functional balance between root and shoot activity, (ii) coordination of leaf photosynthesis. The model was parameterized from plant functional trait measurements of 13 native perennial pasture grass species grown in monocultures at high N availability and low cutting frequency in a field trial. Predicted and measured annual dry-matter yields were highly correlated without bias across species, N supply and cutting frequency treatments in monocultures and in mixtures of six species. Results show the ability of this mechanistic model to simulate without bias nitrogen and disturbance responses of net primary productivity and of plant community structure.

Introduction

How global changes may influence the dynamics of biodiversity across various levels of biological organization is a key challenge to predict future ecosystem functions and services. Ecosystem response to any environmental change is not only driven by the direct effects of abiotic controls but also indirectly affected by changes in the physiology and morphology of individual plants and by the structure of plant populations and communities [START_REF] Suding | Scaling environmental change through the community-level: a trait-based response-and-effect framework for plants[END_REF][START_REF] Klumpp | Using functional traits to predict grassland ecosystem change: a mathematical test of the response-and-effect trait approach[END_REF]. How individual response scale into ecosystem level is sometimes well-known, e.g. photosynthesis scales from leaf to ecosystem [START_REF] Field | Responses of terrestrial ecosystems to the changing atmosphere: a resource-based approach[END_REF]. In contrast, many population (e.g. self-thinning, plasticity) and community processes (e.g. species interactions, species replacement) are not well understood [START_REF] Mcgill | Rebuilding community ecology from functional traits[END_REF][START_REF] Van Wijk | Predicting ecosystem functioning from plant traits: results from a multi-scale ecophysiological modelling approach[END_REF][START_REF] Gross | Linking individual response to biotic interactions with community structure: a traitbased framework[END_REF]. An urgent goal for global environmental change research is thus to better understand the consequences of these complex processes at the population and community levels and how they may impact ecosystem functioning [START_REF] Suding | Scaling environmental change through the community-level: a trait-based response-and-effect framework for plants[END_REF].

Plant functional traits, which are morphological, chemical, physiological and phenological plant characteristics [START_REF] Lavorel | Plant functional classifications: from general groups to specific groups based on response to disturbance[END_REF][START_REF] Violle | Let the concept of trait be functional![END_REF], have been proposed as a useful tool to upscale individual response into ecosystem level [START_REF] Lavorel | Plant functional types: are we getting any closer to the Holy Grail?[END_REF][START_REF] Suding | Scaling environmental change through the community-level: a trait-based response-and-effect framework for plants[END_REF]. Functional traits are linked with individual tolerances to simple abiotic factors [START_REF] Suding | Relationships among species traits: separating levels of response and identifying linkages to abundance[END_REF][START_REF] Louault | Plant traits and functional types in response to reduced disturbance in a semi-natural grassland[END_REF]Gross et al., 2007a) and biotic interactions [START_REF] Soussana | Competition for resources between neighbouring species and patch scale vegetation dynamics in temperate grasslands[END_REF]Gross et al., 2007b[START_REF] Gross | Linking individual response to biotic interactions with community structure: a traitbased framework[END_REF][START_REF] Lavorel | Plant functional types: are we getting any closer to the Holy Grail?[END_REF]. Moreover, functional traits are used to translate individual responses into the community and the ecosystem levels [START_REF] Schymanski | Optimality as a concept to understand and model vegetation at different scales[END_REF][START_REF] Suding | Scaling environmental change through the community-level: a trait-based response-and-effect framework for plants[END_REF][START_REF] Gross | Linking individual response to biotic interactions with community structure: a traitbased framework[END_REF]. At the ecosystem level, the same traits (community mean-trait, [START_REF] Violle | Let the concept of trait be functional![END_REF] determine ecosystem functioning through additive (community effect is determined by mass-ratio processes, [START_REF] Klumpp | Using functional traits to predict grassland ecosystem change: a mathematical test of the response-and-effect trait approach[END_REF] and non-additive (community effect is determined by complementarity processes, [START_REF] Villeger | New multidimensional functional diversity indices for a multifaceted framework in functional ecology[END_REF] effects.

Nevertheless, since almost all empirical studies with plant traits are correlative, the causal linkages between different levels of organizations remain unclear (see review of [START_REF] Lavorel | Plant functional types: are we getting any closer to the Holy Grail?[END_REF][START_REF] Savage | A general multi-trait-based framework for studying the effects of biodiversity on ecosystem functioning[END_REF]. A mechanistic approach could help improve our understanding of the role of plant traits for community dynamics and ecosystem functioning [START_REF] Loreau | Perspectives and challenges[END_REF]. Ideally, this approach would assemble within the same framework 'the explicit inclusion of organismal trade-offs, of environmental constraints, and of the basic mechanisms of interspecific interactions' [START_REF] Tilman | Constraints and tradeoffs: toward a predictive theory of competition and succession[END_REF]. Such requirements imply developing a novel modeling framework by assembling at least three components: (i) biophysical laws that simulate the energy, carbon and water exchanges between vegetation and atmosphere at ecosystem scale (e.g. [START_REF] Wohlfahrt | A model of whole plant gas exchange for herbaceous species from mountain grassland sites differing in land use[END_REF], (ii) carbon, nitrogen and water cycles in the soil-plant-atmosphere continuum (e.g. [START_REF] Thornley | Plant and Crop Modelling: A Mathematical Approach to Plant and Crop Physiology[END_REF], and (iii) plant population dynamics based on resource competition (e.g. [START_REF] Tilman | Plant Strategies and the Dynamics and Structure of Plant Communities[END_REF].

An individual based approach may offer a valuable perspective to build dynamical structure-function-diversity models for several reasons [START_REF] Grimm | A standard protocol for describing individual-based and agent-based models[END_REF]. Across the hierarchy of biological organization levels, selection occurs at the individual level [START_REF] Marks | A holistic tree seedling model for the investigation of functional trait diversity[END_REF]. In addition, individual plants are well defined, measurable and they can be characterized by parameters derived from measured functional traits, that have been recently used for new modeling development in ecology [START_REF] Lehsten | Turnover of plant trait hierarchies in simulated community assembly in response to fertility and disturbance[END_REF][START_REF] Savage | A general multi-trait-based framework for studying the effects of biodiversity on ecosystem functioning[END_REF][START_REF] Van Wijk | Predicting ecosystem functioning from plant traits: results from a multi-scale ecophysiological modelling approach[END_REF]. Moreover, the individual based approach offers the possibility to simulate plastic adjustments of plant form and function [START_REF] Yin | Simulating the partitioning of biomass and nitrogen between roots and shoot in crop and grass plants[END_REF][START_REF] Hoglind | Timothy growth in Scandinavia: combining quantitative information and simulation modelling[END_REF] in response to resource levels mediated by interactions with neighbors.

To date, most dynamic structure-function models based on plant traits concern single trees [START_REF] Le Roux | Carbon-based models of individual tree growth: a critical appraisal[END_REF][START_REF] Allen | Using L-systems for modeling source-sink interactions, architecture and physiology of growing trees: the L-PEACH model[END_REF][START_REF] Marks | A holistic tree seedling model for the investigation of functional trait diversity[END_REF]. Recently, theoretical models have been developed to demonstrate the role of traits, their diversity, their degree of correlation and their plasticity for ecosystem functioning [START_REF] Norberg | Phenotypic diversity and ecosystem functioning in changing environments: a theoretical framework[END_REF][START_REF] Savage | A general multi-trait-based framework for studying the effects of biodiversity on ecosystem functioning[END_REF]. However, these models include a few traits only and do not take into account biophysical laws that constrain mass and energy exchange in plant canopies. Dynamic models of pasture grasses [START_REF] Hoglind | Timothy growth in Scandinavia: combining quantitative information and simulation modelling[END_REF][START_REF] Groot | An object-oriented model of the morphological development and digestibility of perennial ryegrass[END_REF][START_REF] Lafarge | A modelling of the tillering capable of reproducing the fine-scale horizontal heterogeneity of a pure grass sward and its dynamics[END_REF][START_REF] Tomlinson | A functional-structural model for growth of clonal bunchgrasses[END_REF] were developed for monocultures. However, these models cannot be used for mixtures and are sometimes restricted to the simulation of a single growing season [START_REF] Groot | An object-oriented model of the morphological development and digestibility of perennial ryegrass[END_REF]. Finally, non individual-based models simulating mixed grasses [START_REF] Schippers | Competition for light and nitrogen among grassland species: a simulation analysis[END_REF] and grass-clover growth [START_REF] Lazzarotto | Dynamics of grass-clover mixtures: an analysis of the response to management with the productive grassland simulator (Prograss)[END_REF] were previously developed, but they are not based on plant traits framework and they do not include shoot and root morphogenesis.

We have developed from previous works (Soussana et al., 2000a,b) a dynamic structure-function model which simulates average individual plants for each population of a multispecies canopy consisting of perennial C 3 grass species. This model called Gemini (Grassland Ecosystem Model with INdividual centred Interactions) is parameterized from a large number of shoot and root traits in each plant population (Table 1). The aims of this structure-function-diversity model are to understand how plant traits interact with abiotic factors to control plant population dynamics, plant community dynamics and ecosystem functioning. In the present paper, we review the key concepts of this model and discuss its parameterization from plant trait measurements. We provide evidence of the model ability to simulate: (i) within and between species variation in above-ground productivity along environmental gradients for grass monocultures, (ii) species dominance in grass mixtures.

Methods

Experimental data

The experimental site used for model evaluation was established in spring 2002 in an upland area of central France (Theix, 45 • 43 N, 03 • 01 E, 870 m a.s.l.) on a granitic brown soil (Cambisol, FAO). The local climate is semi-continental, with a mean annual temperature of 9 • C and a mean annual rainfall of 760 mm. native perennial C 3 grass species that co-occur in mesic permanent grasslands were studied in monocultures: Alopecurus pratensis, Anthoxanthum odoratum, Arrhenatherum elatius, Dactylis glomerata, Elytrigia repens, Festuca arundinacea, Festuca rubra, Holcus lanatus, Lolium perenne, Phleum pratense, Poa pratensis, Poa trivialis, Trisetum flavescens. A Lolium perenne cultivar (Clerpin) was added as a control. Henceforth, in the text, species are referred to by their species name. In addition, three mixtures of six species, which were drawn among the 13 species, were also studied. The mixtures were: (i) D. glomerata, F. arundinacea, F. rubra, L. perenne, P. pratensis and C. cristatus; (ii) A. pratensis, A. odoratum, A. elatius, E. repens, H. lanatus and T. flavescens; (iii) A. elatius, D. glomerata, E. repens, F. arundinacea, F. rubra and H. lanatus. C. cristatus could not be studied in monoculture as it suffered freezing damage at the time of seedling implantation in monoculture during the winter. This species represented very low abundance (<2%) within mixture and was not considered in this study.

The experimental design has 3 complete randomized blocks each crossing two factors: cutting frequency (3 and 6 cuts yr -1 , C-and C+, respectively) and N fertilizer supply (120 and 360 kgN ha -1 yr -1 , N-and N+, respectively). Phosphorus and potassium were supplied in spring at non-limiting rates for growth. When soil water content (SWC) was below 10%, all plots were irrigated (see [START_REF] Pontes | Leaf traits affect the above-ground productivity and quality of pasture grasses[END_REF] for full details). Plant functional traits of monocultures (Table 1) were measured in 2003, and 2006 in the high N supply (N+) and low disturbance (C-) treatment [START_REF] Maire | Trade-off between root nitrogen acquisition and shoot nitrogen utilization across 13 cooccurring pasture grass species[END_REF][START_REF] Pontes | The role of plant traits and their plasticity in the response of pasture grasses to nutrients and cutting frequency[END_REF]. This treatment provided non-limiting conditions for morphogenetic development and above-ground productivity; and was used for model parameterization. Other monoculture treatments (C+N+, C-N-, C+N-) as well as all mixtures were used for an independent evaluation of the model.

Model purpose

The model is described following the ODD (Overview, Design concepts and Details) standard protocol proposed by [START_REF] Grimm | A standard protocol for describing individual-based and agent-based models[END_REF] for individual-based and agent-based models. It should be noted that Gemini is an individual-centred model, rather than being individual-based, since it simulates average individuals within each plant (or animal) population. A detailed list of all 132 equations, as well as the 187 variables and the 100 default parameter values and their units is available (at www1.clermont.inra.fr/urep/modeles/gemini.htm) and will be send on request.

The main purpose of Gemini is to understand the dynamics and plasticity of plant species within a community and the role of traits and their plasticity for ecosystem functioning. The model considers climatic (short-wave radiation, temperature and precipitation) and atmospheric (CO 2 concentration) abiotic drivers. Management conditions concern both disturbance (by cutting and grazing) and fertilization (inorganic and organic N supply). The model was built Plant traits used for the calibration of the Gemini model. Abbreviations, units and references are given for each trait ([1] [START_REF] Cornelissen | A handbook of protocols for standardised and easy measurement of plant functional traits worldwide[END_REF][2] [START_REF] Kazakou | Components of nutrient residence time and the leaf economics spectrum in species from Mediterranean old-fields differing in successional status[END_REF][3] Pontes et al., 2007;[4] [START_REF] Pontes | The role of plant traits and their plasticity in the response of pasture grasses to nutrients and cutting frequency[END_REF][5] [START_REF] Craine | Functional traits, productivity and effects on nitrogen cycling of 33 grassland species[END_REF]. Each trait was used to calibrate one or a few parameters by species. Parameter definitions are provided in

mg g -1 DM h -1 [1,5] Su0 Physiology Leaf and root chemical composition trait Leaf N content LNC mg g -1 [1-5] fns Chemical composition Leaf C content LCC mg g -1 [5] fcs Chemical composition Root N content RNC mg g -1 [5] fnr Chemical composition Root C content RCC mg g -1 [5] fcr Chemical composition
with a modular architecture, which permits to include or not different biotic agents (plant species, soil microbial decomposers, and herbivores at grazing) as well as different environment and management modules (soil, vegetation, fertilization and cutting).

Gemini can simulate a potentially unlimited number of plant species (or plant populations from the same species) from currently two plant functional types (perennial grasses and legumes). The model focuses on the acquisition and the utilization of two major resources (light and nitrogen) by plants and their consequences for the carbon and nitrogen cycles.

State variables and scales

Gemini consists of vegetation, soil and herbivore sub-models, coupled with environment and management sub-models (Fig. 1). The vegetation sub-model, named Canopt (Soussana et al., 2000a,b) is an individual-centred model of a multi-species stand comprising clonal grasses and/or legumes and forming a multi-layer plant canopy. Each clonal plant population is described as a collection of identical axes (e.g. tillers for grasses). Moreover, all plant species are assumed to be perfectly mixed in the horizontal plane. Plant population demography is calculated from the vegetative multiplication and mortality of axes.

The vegetation sub-model consists of four modules: (i) a biochemical module, which simulates the C and N balance and the partitioning of growth among shoot structures (W S ), leaf proteins (W P ) and roots (W R ) of a collection of identical axes for each plant population. The corresponding state variables are the number of axes by population (D), the masses of three structural compartments, of two C and N substrate compartments and of four C and N reserve compartments; (ii) a shoot morphogenesis module, which computes the demography and size of leaves (two state variables, length and mass per leaf); (iii) a root morphogenesis module, which computes the demography and size of roots (two state variables, length and mass, per root); (iv) a competition module which calculates short-wave radiation and inorganic N partitioning among mixed plant species.

The environment sub-model calculates the microclimate within the canopy and the inorganic N balance of the soil (or of the substrate when the vegetation model is not coupled to the soil model).

The management sub-model schedules events caused by grassland management (cutting dates, grazing periods, N fertilizer applications).

Process overview and scheduling

The carbon balance is based on [START_REF] Farquhar | A biochemical model of photosynthetic CO2 assimilation in leaves of C3 species[END_REF] equations for leaf photosynthesis which were updated according to [START_REF] Maire | From functional traits of grasses to the functioning of grassland ecosystem: A. Mechanistic modelling approach[END_REF]. A linear relationship of Vc max (the maximal carboxylation activity of Rubisco) and J max (the electron transport capacity) with the area based leaf protein concentration [START_REF] Field | Allocating leaf nitrogen for the maximization of carbon gain: leaf age as a control on the allocation program[END_REF][START_REF] Nijs | Leaf nitrogen content as a predictor of photosynthetic capacity in ambient and global change conditions[END_REF] is assumed. The vertical profile of leaf proteins is simulated according to [START_REF] Hirose | Maximizing daily canopy photosynthesis with respect to the leaf nitrogen allocation pattern in the canopy[END_REF]. Respiration is divided into growth and maintenance components. Leaf respiration varies with leaf protein content.

Nitrogen acquisition is based on soil diffusion processes from soil to roots surface [START_REF] Barber | Soil Nutrient Bioavailability. A Mechanistic Approach[END_REF]. The model separates metabolically active roots and roots which have entered a first stage of senescence (see morphogenesis) and have lost their uptake capacity.

The relative growth rate of the three structural compartments (W P , W S , W R ) is calculated through a bi-substrate C-N growth equation using three partitioning variables, one for each compartment [START_REF] Johnson | A model of shoot:root partitioning with optimal growth[END_REF]. The partitioning sub-model specifies a target root/shoot ratio, according to the functional balance hypothesis [START_REF] Brouwer | Nutritive influences on the distribution of dry matter in the plant[END_REF][START_REF] Davidson | Effect of root/leaf temperature differentials on root/shoot ratios in some pasture grasses and clover[END_REF][START_REF] Hilbert | A model allocating growth among leaf proteins, shoot structure, and root biomass to produce balanced activity[END_REF]

Canopt

Plant can Shoot

Leaf Root System Root

1 1 1 1 1 * 0..1 * * 0..1

Leaf of a fun ctional group

Shoot of a functional gr oup Functional g roup

A B

A is compos ed by B B is a specialisation of A A B 1991), and a target protein/structure ratio for leaves, according to the coordination theory of leaf photosynthesis [START_REF] Chen | Coordination theory of leaf nitrogen distribution in a canopy[END_REF][START_REF] Maire | From functional traits of grasses to the functioning of grassland ecosystem: A. Mechanistic modelling approach[END_REF]. The C and N substrates are immediately available both for shoot and root developments, although C and N reserves are assumed to follow first order kinetics. The potential leaf extension rate is calculated according to [START_REF] Van Loo | Tillering, leaf expansion and growth of plants of two cultivars of perennial ryegrass grown using hydroponics at two water potentials[END_REF]. This rate is assumed to increase exponentially with temperature [START_REF] Lemaire | Leaf tissue turnover and efficiency of herbage utilization[END_REF]. The duration of leaf extension is proportional to the phyllochron (leaf emission per unit thermal time) and to the number of leaves growing simultaneously on the same axis. The potential leaf area and leaf weight are calculated from leaf length through simple allometric relationships. Whenever the supply of assimilates to the shoots does not allow this potential [START_REF] Poorter | Construction costs and payback time of biomass: a whole plant perspective[END_REF], the leaf extension rate is reduced and the phyllochron is increased. Hence, the model considers a source limitation of shoot growth when assimilates (C and N) supply to the shoots is less than their potential growth rate. The leaf lifespan is calculated from the maximum number of green leaves per axis [START_REF] Lemaire | Leaf tissue turnover and efficiency of herbage utilization[END_REF]. Roots decay according to a first order kinetics, based on thermal time. Fixed fractions of the shoot structures, leaf proteins and roots are recycled to the C and N substrate compartments during senescence.

GEMINI CANO PT

Plants are assumed to consist of a collection of identical and anatomically connected axes (e.g. tillers for grasses). A maximum site-filling is assumed [START_REF] Davies | Leaf tissue remaining after cutting and regrowth in perennial ryegrass[END_REF][START_REF] Neuteboom | Tillering potential and relationship between leaf and tiller production in perennial Ryegrass[END_REF][START_REF] Van Oijen | Process based modelling of timothy regrowth[END_REF] to calculate the potential branching (or tillering) rate per unit thermal time. Branching is further constrained by light transmission at ground level [START_REF] Simon | Tillering and leaf area index in grasses in vegetative phase[END_REF][START_REF] Bouman | Description of the growth model LINGRA as implemented in CGMS[END_REF] and by the amounts per plant of C and N substrates which are available for the growth of new branches (or tillers). According to the degree of physiological integration in the clonal plant, it is considered that only a fraction of C and N substrates are made available for branching. The baseline mortality rate of axes is calculated from a first order kinetics based on thermal time. This baseline mortality is increased whenever the relative growth rate becomes negative. In this case, mortality rate increases when labile reserves have a low level. During axis senescence, fixed fractions of shoot structures, leaf proteins and roots are recycled to the C and N substrate compartments. At each cutting date, it is assumed that all the plant material above cutting height is removed from the patch. The reserve compartments are assumed to be located below the sward defoliation height.

In order to compute a radiative balance, the canopy is divided into fixed height horizontal layers and the light absorbed by leaves from each species in each layer is calculated assuming: (i) a full diffuse radiation, (ii) a perfect mixing of leaves in each horizontal layer. The daily photosynthetically active radiative balance of each population in the mixture is calculated following [START_REF] Sinoquet | Comparison of models for daily light partitioning in multispecies canopies[END_REF] using Kubelka-Munk equations for horizontally homogeneous multispecies canopies. Canopy photosynthesis is then calculated for each population by cumulating the leaf photosynthesis of each canopy layer.

Competition for inorganic N between mixed plant populations is calculated from the total demand for N uptake (estimated from the uptake capacity and from the cumulated root length of each population). A diffusion approach [START_REF] Sheehy | A nitrogen-led model of grass growth[END_REF] is used to calculate the daily amount of N supplied by the soil to the roots. Whenever N supply per unit root area is lower than the uptake capacity, the actual uptake rate is calculated by considering that each plant population absorbs the same fraction of its potential demand per unit root area.

Software concepts

The first step of a simulation (Fig. 2) is the construction of the bio-physical model (Part 1 simulation start). The hierarchical structure of the object-oriented code allows for a number of possibilities, including one to n plant species; describing explicitly or not, the root and shoot morphogenesis; including or not one to n soil layers and an animal object. Then, for each simulation day, the main loop (Part 2) goes through three successive steps: (i) update light and N availability in each horizontal soil and canopy layer; (ii) allocate resources among plant species; (iii) update all objects at the end of simulation day for each species. At the end of the daily loop, the simulation goes one step forward if at least one plant species has survived. The loss of a plant species occurs either when population axis density becomes nil, or when the simulated mean axis has no more substrates.

The model runs on a robust, flexible and portable platform. It allows coupling different numerical sub-models composed by sets of ordinary differential equations [START_REF] Cellier | Continuous System Simulation[END_REF]. This software platform, UNIF (Unified Numerical Integration Framework), is based on a fully object-oriented structure. It has been developed with the C ++ language and uses the classical relationships of the object-oriented paradigm (mainly inheritance and composition). These objects are referenced in runtime modifiable lists. To check for numerical quality, mass balance (C and N) is continuously evaluated by independent calculations. Simulations of a plant population were stopped whenever the substrate C concentration increased above an unrealistic value (greater than the structural mass). Moreover, UNIF allows pausing simulation and saving all state variables of user defined objects (e.g. plant populations). This allows for spin-up runs of plant species in monocultures, prior to the construction of a more diverse community. This restart procedure allowed initial plant populations to be close to steadystate.

Design concepts

Coordination of C and N fluxes

Coordination of leaf and canopy photosynthesis. The formulation and parameterization used by [START_REF] Maire | From functional traits of grasses to the functioning of grassland ecosystem: A. Mechanistic modelling approach[END_REF] for the biochemical of leaf photosynthesis by [START_REF] Farquhar | A biochemical model of photosynthetic CO2 assimilation in leaves of C3 species[END_REF] is used by the model. Assimilation rate is usually limited either by the carboxylase activity of the Rubisco or by the electron flux through the chloroplast photosystems. The maximal rates of these two processes (Vc max and J max , respectively) scale linearly with leaf protein content [START_REF] Wohlfahrt | A nitrogen sensitive model of leaf CO2 and H2O vapour gas exchange: application to 13 key species from differently managed mountain grassland ecosystems[END_REF][START_REF] Wohlfahrt | Inter-specific variation of the biochemical limitation to photosynthesis and related leaf traits of 30 species from mountain grassland ecosystems under different land use[END_REF], which led to the coordination theory [START_REF] Chen | Coordination theory of leaf nitrogen distribution in a canopy[END_REF]. According to this theory, under given environmental conditions, there is a single leaf N content which results in a co-limitation of leaf photosynthesis by Rubisco activity and by RuBP regeneration. This area-based leaf photosynthetic protein content (N pac , gN m -2 ) maintains a balance between the two processes potentially limiting photosynthesis. Below N pac , the photosynthesis will be limited by the Rubisco activity and therefore by the amount of leaf proteins. Beyond N pac , the marginal gain of photosynthesis per unit of leaf proteins is weak. In the Gemini model, the coordinated leaf N content is a target concentration for the allocation of assimilates between leaf structure and leaf photosynthetic proteins. This allows calculating a single leaf photosynthetic N content per unit area (N pac ), which ensures photosynthesis co-limitation by light and by CO 2 capture and use (Soussana et al., 2000a).

C and N substrates allocated to leaf proteins are distributed among horizontal canopy layers in proportion to leaf area and to the coordinated leaf protein content per unit leaf area (N pac ). Along the vertical canopy profile, N pac declines with transmitted PPFD when all other variables are equal. The coordination theory matches the predicted linear relationship between photosynthetic proteins and transmitted PPFD [START_REF] Chen | Coordination theory of leaf nitrogen distribution in a canopy[END_REF] which has often been reported [START_REF] Lötscher | Vertical leaf nitrogen distribution in relation to nitrogen status in grassland plants[END_REF]. Air temperature (T a ), relative humidity (h s ) and CO 2 concentration (C a ) also vary along the vertical canopy profile. Higher h s and lower T a at depth would reduce N pac for a given PPFD, while a lower C a would increase it. However, variations in T a , h s and C a along the vertical canopy profile are currently neglected by the model.

Coordination of C-N fluxes at plant scale. In higher plants, autotrophy is obtained by combining light capture by shoots and nutrients and water capture by roots. For an individual plant, the most efficient partitioning of growth would ensure that no resource is in excess, which implies that no single resource becomes limiting [START_REF] Chapin | Integrated responses of plants to stress. A centralized system of physiological responses[END_REF]. A first version of this multiple resource capture hypothesis is the functional balance hypothesis which assumes that partitioning of growth between shoots and roots tends to balance shoot photosynthesis and root N acquisition to ensure C to N homeostasis. This hypothesis was shown to account for a large number of empirical observations (see the review by [START_REF] Wilson | A review of evidence on the control of shoot:root ratio, in relation to models[END_REF].

In an original way, Gemini combines the functional balance and photosynthetic coordination theories to determine partitioning rules for potential assimilate supply to the three structural compartments. Following [START_REF] Hilbert | A model allocating growth among leaf proteins, shoot structure, and root biomass to produce balanced activity[END_REF], P is the partitioning ratio between shoots structures and roots. Q is the partitioning ratio between shoot structures and leaf proteins. Therefore, the functional balance hypothesis is the main determinant of the balance between shoot structures and roots, while the coordination hypothesis of leaf photosynthesis determines the potential partitioning among leaf proteins and shoot structures. By varying the partitioning ratios between shoots and roots, and between shoot structures and leaf proteins, the model converges towards a co-limitation of plant growth by the captures of light, CO 2 and N.

Regulation of root N capture. The model considers that both morphological and physiological plasticity regulate N capture. First, a hierarchy of root orders with different diameters is considered by the model and differential root growth and branching leads to increased root area per unit mass in case of N deprivation. Second, the uptake capacity may become down-regulated whenever N substrate concentration within the plant increases above a fixed threshold [START_REF] Imsande | N demand and the regulation of nitrate uptake[END_REF][START_REF] Soussana | A simple model of feedback regulation for nitrate uptake and N2 fixation in contrasting phenotypes of white clover[END_REF].

Internal C and N recycling

Storage and remobilization. Pasture grasses accumulate reserve carbohydrates such as fructans and other glucidic polymers mostly within the vacuole. These carbon reserves tend to escape defoliation since they are concentrated within the sheath [START_REF] Millard | The accumulation and storage of nitrogen by herbaceous plants[END_REF][START_REF] Jeuffroy | Integrated physiological and agronomic modelling of N capture and use within the plant[END_REF].

Gemini considers labile and slow reserve C and N pools. Slow reserve pools are especially important in species developing specialized organs such as rhizomes (e.g. Elytrigia repens) and tap roots (e.g. Festuca paniculata). Both labile and slow reserves are modeled using a buffer formalism that can be tuned to increase remobilization intensity at the expense of the duration of the remobilization period.

Resorption at senescence. Gemini follows the efficiency hypothesis [START_REF] Franklin | Leaf senescence and resorption as mechanisms of maximizing photosynthetic production during canopy development at N limitation[END_REF][START_REF] Kazakou | Components of nutrient residence time and the leaf economics spectrum in species from Mediterranean old-fields differing in successional status[END_REF], by considering distinct efficiencies of resorption for leaf proteins, for carbon and nitrogen both in roots and in shoot structures. Substrate and reserve C and N undergo full resorption during senescence. Hence, C and N in labile pools (substrates, reserves, and proteins) are assumed to be mostly recycled during senescence, while cell walls from leaf structures are mostly returned to the soil.

Morphogenesis model

Plant morphology and architecture. Clonal plants can be seen as colonies of population units, i.e. tillers in grasses [START_REF] Gillet | Les graminées fourragères. Description, fonctionnement[END_REF]. During the vegetative growth of grasses, a single meristem forms all shoots of a population unit, which can therefore best be described as an axis. This axis forms several primary root meristems, usually two per leaf emitted [START_REF] Matthew | Determination of tiller and root appearance in perennial ryegrass (Lolium perenne) swards by observation of the tiller axis, and potential application in mechanistic modelling[END_REF]. These meristems develop primary roots that may develop secondary meristems leading to the outgrowth of secondary roots. The maximum root branching order is usually less than three in pasture grasses.

This simplified description of plant populations and their morphology is used by the Gemini model which effectively combines three hierarchical levels, the individual root or leaf, the axis and the plant population. Based upon the morphological description at the axis level, the model calculates 2D architecture by defining the relative positions of the bases and tips of each root and leaf. The relative positions of leaf bases are all located at the soil surface. The azimuth of roots and of leaves is not calculated; therefore the model provides only a 2D description and not a full 3D architecture. Roots are assumed to be linear segments, which length is however greater or equal to the distance between tip and basis, since root tortuosity parameter value can be greater than one. The shape of individual grass leaves is assumed to follow a rectangular hyperbola, which form varies according to the canopy height and tiller density (Soussana et al., 2000b).

This axis-centered 2D architectural description is then scaled up to the population level by assuming that all plant axes are identical, thereby generating a 1D distribution of leaf and of root area in the horizontal canopy and soil layers which are used by the competition sub-model. This simplified approach of plant morphology does not account for stems and for internodes which may position leaf bases in upper canopy layers. Therefore, the morphological sub-model is restricted to vegetative grasses.

Coordination of potential morphogenesis. A general theory of resource distribution through hierarchical branching networks has been applied to the case of vascular plants [START_REF] West | A general model for the structure and allometry of plant vascular systems[END_REF], showing the role of the topology of the branching network for the assimilate demand and for the body size [START_REF] Enquist | Global allocation rules for patterns of biomass partitioning in seed plants[END_REF]. In Gemini, morphogenesis generates a simplified branching network of shoot and root parts through six categories of morphogenetic processes: emission, elongation, vegetative multiplication, branching, senescence and axis death. Each of these processes has a potential temperature dependent rate. Leaf elongation and vegetative multiplication are, in addition, controlled by light-mediated signals (photoperiod and light extinction at ground level, respectively). The number of growing leaves and roots is calculated each day, based on the thermal time between emission and maturity [START_REF] Lemaire | Les flux de tissus foliaires au sein des peuplements prairiaux. Eléments pour une conduite raisonnée du pâturage[END_REF]. For each growing plant part, potential morphogenesis generates each day an additional volume, through an allometric relationship with the length increment generated by elongation. Hence, the total volume increment generated by shoots and roots growth can be calculated on a daily basis per plant axis. Some morphogenetic processes are coordinated. The thermal time between two successive emission events of leaves or roots (i.e. the phyllochron for leaves) is used as a biological clock controlling the ontogeny of the plant. First, the duration of leaf elongation varies with phyllochron [START_REF] Lemaire | Les flux de tissus foliaires au sein des peuplements prairiaux. Eléments pour une conduite raisonnée du pâturage[END_REF], since the number of growing leaves per axis is assumed to be constant (between 1 and 2 for grasses, [START_REF] Gillet | Les graminées fourragères. Description, fonctionnement[END_REF]. Second, leaf senescence starts after a thermal time which is a multiple of the phyllochron since the number of mature leaves per axis is assumed to be constant. Third, vegetative multiplication is also constrained by phyllochron, since the potential site-filling rule [START_REF] Neuteboom | Tillering potential and relationship between leaf and tiller production in perennial Ryegrass[END_REF] states that the number of axillary buds varies with leaf number and, hence, with phyllochron.

While roots and leaves share common properties (in terms of emission and elongation), roots also have specific properties. Root hierarchy constrains the emission rate, the diameter, length and longevity of roots. Roots of higher order are emitted more frequently, are thinner, shorter and short-lived compared to roots of lower branching order [START_REF] Pagès | Root Typ: a generic model to depict and analyse the root system architecture[END_REF].

Root branching only occurs for roots with an order lower than the maximum (e.g. when the maximum branching order is 3, only primary and secondary roots can branch). The model implicitly assumes that a pair of primary roots is emitted for each phyllochron. However, the primary roots turnover is not explicitly simulated on the average axis, since the primary roots death corresponds to axis death, which is simulated through the axis mortality process. Rather, axis death is calculated at the population level, through a baseline death rate which is increased whenever relative growth rate becomes negative and when reserves are low.

Leaf and root (of order greater than one) senescence affects the area of these plant parts. With leaves, senescence starts from the leaf tip and progresses towards the leaf basis. Dead leaves are assumed to be detached from the plant axis. Root senescence is also a gradual process, with three successive stages: (i) transition from metabolically active roots to roots which have lost active metabolism but can still transport water, (ii) transition from 'passive' roots to dead but still attached roots, (iii) dead and detached roots. Leaf and root senescence tend to balance emission, thereby leading to a steady-state architecture which converges towards a constant number of growing, mature and senescing leaves and of active, passive and dead primary and higher order roots. Therefore, potential rates of the morphogenetic processes constrain the morphology and architecture of the plant populations, leading to coordination between leaf and root numbers.

Model parameterization

Out of a total of 100 parameters, 65 have been considered as constant across the 13 perennial grass species investigated with a value taken from the literature. The remaining parameters (Table 2) were calculated from plant trait values measured in the field experiment comparing 13 grass species grown in monocultures in the C-N+ management treatment. In this treatment, neither nutrients nor water resources were limiting above-ground grass growth [START_REF] Pontes | Leaf traits affect the above-ground productivity and quality of pasture grasses[END_REF][START_REF] Maire | Trade-off between root nitrogen acquisition and shoot nitrogen utilization across 13 cooccurring pasture grass species[END_REF].

Shoot morphology parameters were derived from single trait measurements (for L leaf0 , maximal leaf length, F sheath , fraction of sheath in total leaf length; ng leaf , nm leaf , growing and mature leaf number per tiller, ph 0 , phyllochron and LDMC, leaf lamina dry-matter per unit fresh-matter), or from combined trait measurements (for leaf allometry parameters a AL , b AL , a WL , b WL were fitted from relationships between lamina area or volume and leaf length). The last two remaining shoot parameters (C, which controls leaf lamina angle for a given canopy height and density, and b LER , power coefficient of temperature response of leaf elongation) were estimated from measurements of leaf angles and from observations of earliness of vegetative spring growth, respectively.

Root morphology parameters were derived from measurements on root ingrowth cores filled with vermiculite and supplied with slow release fertilizers (for full details see [START_REF] Maire | Trade-off between root nitrogen acquisition and shoot nitrogen utilization across 13 cooccurring pasture grass species[END_REF]. In the same way, parameter values were derived from single trait measurement (rr root1 , radius of primary roots, RDMC, root drymatter content, Order max , maximum root branching order), or from combined traits measurements (Orl, branching order dependency, which were obtained from radius measurements of root orders). L rootmax1 , the potential length of primary roots was estimated from [START_REF] Kutschera | Wurzelatlas mitteleuropäischer ackerunkräuter und kulturpflanzen[END_REF] root atlas and from the correlation between maximal plant height and rooting depth [START_REF] Craine | Functional traits, productivity and effects on nitrogen cycling of 33 grassland species[END_REF].

Two parameters (LL plast , the branching order dependency of potential root length and Tr, the thermal time interval between two successive root emission events) were optimized by maximizing axis biomass (W G ). This first optimization was done on C-N+ management treatment using constant axis density for each species. The two population demography parameters (intcl, the clonal integration and Tsen 0 , the lifespan of an axis) were then optimized by fitting simulated to measured tiller density (D) per unit ground area. This second optimization was done after the root parameters optimization.

Root chemical composition parameters were derived from measured C and N contents in below-ground plant parts. C fraction in shoot structures was derived from measurements in mature green leaves. The N fraction of shoot structures was estimated from the measured N contents of senescent leaves [START_REF] Maire | Trade-off between root nitrogen acquisition and shoot nitrogen utilization across 13 cooccurring pasture grass species[END_REF].

Finally, species specific physiological parameters measurements were derived from measurements of (i) N (R P ) and dry-matter (R C , R N ) resorption during leaf senescence [START_REF] Maire | Trade-off between root nitrogen acquisition and shoot nitrogen utilization across 13 cooccurring pasture grass species[END_REF], (ii) root N uptake potentials (Su 0 ) in root ingrowth cores [START_REF] Maire | Trade-off between root nitrogen acquisition and shoot nitrogen utilization across 13 cooccurring pasture grass species[END_REF]. The potential size of the slow C reserve pool (W rrCmax ) was estimated for each species as the minimum size stabilizing the C substrate pool during winter. The potential size of the slow N reserve pool was estimated at one tenth of that of the slow C reserve pool.

Model evaluation

The model was then evaluated by comparing simulated to observed above-ground DM productivity in 2004, two years after the start of experiment, with monocultures of the 13 grass species at low cutting frequency (C-) and at low N under high cutting frequency (C+,N-), that is for experimental treatments which were not used for model parameterization. Finally, the model was evaluated by comparing simulated and observed species relative abundance in three mixtures of six out of the 13 species (see [START_REF] Pontes | Impacts of species interactions on grass community productivity under contrasting management regimes[END_REF].

In monocultures, the restart procedure was used to initialize all simulations. In mixtures, an initial total density of 1000 axes m -2 (equally distributed among the six species) was assumed. For comparison with field data, simulations were run with daily climatic data recorded at the experimental site in 2003 and 2004 (temperature, PPFD, precipitation and air humidity). The same timing (day of year) was applied for the fertilization and cutting events in the model and in the experimental treatments, as well as the same amounts of fertilizer N supply. Since there was no soil assumed, a daily net soil mineralization was assumed at 0.11 gN m -2 d -1 . The dynamics of N fertilizer release in the soil were simulated, assuming a pulse and decline over one month.

The model was evaluated by using four complementary indices [START_REF] Willmott | Some comments on the evaluation of model performance[END_REF][START_REF] Stöckle | Evaluation of estimated weather data for calculating Penman-Monteith reference crop evapotranspiration[END_REF] which compare simulated and measured above-ground DM yields of species monocultures and species abundance in mixtures. The index of agreement (d), which is ranged between 0 and 1 (1 indicating perfect agreement), gives the standard error variability (Eq. (I)). A d value above 0.9 was considered as acceptable error. The relative root mean square error (RE), which is bounded by 0 indicating perfect agreement between simulated and measured values, gives the general standard deviation from the measured mean by placing, however, relatively more weight on high quadratic error values (Eq. (II)). A RE value below 0.2 was considered as acceptable error. Since RE does not specify the error source, a ratio between systematic (RE S ) and unsystematic (RE U ) errors was also calculated according to (Eq. (III)).

d = 1 - n i=1 (E i -M i ) 2 n i=1 (|E i -M| + |M i -M|) 2 (I) RE = n i=1 (E i -M i ) 2 n 0.5 • 1 M (II) RE S RE U = n i=1 ( Êi -M i ) 2 n i=1 (E i -Êi ) 2 with Êi = b • M i + a (III)
where E i and M i are the simulated and measured values of one variable of species i, Mi is the average of M i and Êi is an estimate of E i , which is derived from the linear regression between E i and M i .

The slope of reduced major axis regressions (RMA slope) between estimated and measured values is also an index of model under-or-over-prediction trends. Finally, to test the model ability to simulate cross-species differences, simulated and measured Spearman's ranks of species were compared. Rank tests were performed

Results

In the C-N+ treatment of monocultures, four parameters (i.e. two from root morphology, OrL and Tr, and two from axis demography, Tsen 0 and Intcl) had to be optimized from measured annual above-ground DM yield (DMY) and tiller density, respectively. Moreover, most values of species specific parameter were derived from plant trait measurements in this field treatment (Table 1). Therefore, the good agreement (d = 0.96, RE = 0.16, Table 3A) between simulated and measured DMY means in the C-N+ treatment merely indicates that the optimization procedure was successful but cannot be used to assess the model accuracy.

The three other management treatments were then simulated by Gemini using exactly the same set of parameter values and modifying only the cutting frequency and N supply factors. Despite a lower agreement between observations and simulations for d and RE in the C+N+ and C-N-treatments respectively (indicating a higher variability of model error and a higher quadratic error, respectively), simulated and measured annual DMY of the three independent treatments were overall in good agreement (d ≥ 0.9, except for C+N+ treatment; RE ≤ 0.2 except for the C-N-treatments) with no systematic bias (RE S /RE U < 1) and no under-or-over estimation (RMA slope never different from 1) (Fig. 3A, Table 3A). Moreover, plant species ranks were significantly correlated for simulated and measured DM yields of the grass monocultures (Fig. 3A, insert and Table 3A).

Within species mixtures, simulated and measured species abundances agreed (d ≥ 0.9) and their ranks were significantly correlated. In addition, there was neither a systematic bias (RE S /RE U < 1) nor an under-or-over estimation (RMA slope not significantly different from 1) (Fig. 3B, Table 3B). However, a low overall RE value (0.51) indicated a few large quadratic errors between simulated and measured abundances (RE > 0.2). While the relative abundances of dominant (Fig. 3B, insert; D. glomerata, F. arundinacea, H. lanatus were above 20%) and subordinate (A. elatius, P. pratensis, A. pratensis were above 10%) species were correctly simulated, the model showed a tendency to over-estimate the relative abundance of rare (relative abundance below 10%) species.

Discussion

By combining an individual centered approach coordinating C and N metabolisms, shoot and root morphogenesis and tiller demography, we have developed a mechanistic model capable to link plant functional traits, plant populations and community dynamics, and ecosystem fluxes (e.g. net primary productivity). Despite its apparent complexity, the Gemini model is mostly based on sets of widely published and well established mechanistic equations which are derived from first order principles. Moreover, the model brings a high degree of realism, through an explicit dynamic description of plant form and function, and most, if not all, output variables can be compared with ecophysiological, ecological and agronomical measurements both in the field and in controlled environment studies.

While the model requires a large number of parameters for each plant population, we have shown that it can be readily parameterized from root and shoot trait values measured in field experiments. The development of international trait databases (e.g. TRY, [START_REF] Kattge | TRY -a global database of plant traits[END_REF] will expand the list of traits and of species that can be used for model parameterization.

After parameterization under optimal environmental conditions, the model was able to simulate accurately both inter and intra-specific variation in response to nitrogen and disturbance of pasture grass dry-matter yields (Table 3A, Fig. 3A). Indeed, for the three independent management treatments (C-N-, C+N-, C+N+) which were not used for model parameterization: (i) simulated and observed species ranks for DM yields in monocultures were significantly correlated without bias and (ii) acceptable model errors were shown by two statistical indices (d and RE).

The dual description (morphology and metabolism) of plant populations provided by the model is unified considering a minimum law. Both ontogeny (morphogenetic processes) and functional balance (physiological processes) may limit the actual growth rate of roots and shoots. Plant populations are therefore driven through a dually constrained manifold that keeps their form and function within plastic limits. Supply limitation occurs whenever assimilates supply is less than assimilates demand. Demand limitation occurs whenever assimilate demand is less than assimilate supply, leading to excess assimilates being allocated to the C and N substrate pools.

Resorption during leaf and root senescence leads to C and N fluxes towards the substrate pools. Therefore, part of the assimilate supply for tissue growth depends upon tissue senescence. In the same way, axis death provides through resorption additional substrates, which may be used for vegetative multiplication, that is the emission and outgrowth of new axes. Thus, ontogeny not only requires, but also provides assimilates.

In a dynamic loop, morphogenesis leads to changes in the length and area of root and shoot captors, which affect resource capture and availability and hence assimilate supply for morphogenesis. At the population level, vegetative multiplication of axes leads to lower resource availability for each individual axis, to lower substrate levels and, hence, to increasing limitations on the formation of new axes.

Vegetative multiplication implies a transition from heterotrophy to autotrophy of the newly formed axes. Autotrophy is assumed at the stage when the net carbon balance of a newly formed axis (i.e. photosynthesis minus maintenance respiration) becomes positive, as a result of leaf extension into unshaded canopy layers. Thus, autotrophy is gained earlier in open compared to dense stands. On average, autotrophy accounts for two thirds (data not shown) of the Predicting directional changes in grassland community structure in response to environmental factors has proved to be difficult [START_REF] Suding | Relationships among species traits: separating levels of response and identifying linkages to abundance[END_REF]. Gemini simulated adequately the plant community structure within mixtures of six species (Table 3B, Fig. 3B), thereby capturing directional changes in community structure in response to nutrients (N supply) and disturbance (cutting frequency) drivers. Despite this overall agreement, relatively high quadratic errors were found in some instances, indicating that there is scope for further improvements of the model. Given the trade-off reported for the same grass species [START_REF] Maire | Trade-off between root nitrogen acquisition and shoot nitrogen utilization across 13 cooccurring pasture grass species[END_REF] between nitrate and ammonium root uptake capacities, simulating the uptake and soil dynamics of these two mineral N forms could help refining the model and reducing prediction errors of grass species relative abundance (a significant relationship was found between the residuals of the observed/predicted abundance curve and the root uptake capacities of these species, data not shown).

We do not know of another mechanistic grassland model that can predict nitrogen and disturbance effects on relative abundance in mixtures including more than three species (but for three species mixtures, see [START_REF] Van Wijk | Predicting ecosystem functioning from plant traits: results from a multi-scale ecophysiological modelling approach[END_REF]. In recent years there have been attempts to predict changes in vegetation structure caused by global change, with dynamic global vegetation models (DGVMs, [START_REF] Cramer | Global response of terrestrial ecosystem structure and function to CO2 and climate change: results from six dynamic global vegetation models[END_REF]. These models consider several plant functional types (such as evergreen and deciduous trees, C 3 and C 4 grasses) that co-occur within a given grid box, and some models (Hybrid, Friend et al., 1997;LPJ, Sitch et al., 2003;Orchidee, Krinner et al., 2005) simulate competition between these types. The Gemini model uses plant traits rather than plant functional types to simulate both directional changes in species composition and their effects on ecosystem functioning. This approach could be extended within the framework of future DGVMs [START_REF] Van Wijk | Predicting ecosystem functioning from plant traits: results from a multi-scale ecophysiological modelling approach[END_REF].

Conclusion

The Gemini model provides for the first time a trait-based mechanistic framework for modeling the role of grass diversity and plasticity on community dynamics and ecosystem functioning. The aim of this paper was to first present the model and evaluate its behavior. The first results presented here show that the model can predict productivity patterns both within and across species, and simulate the main directional changes in plant community structure in response to environmental factors. In a second paper (Maire et al., in press), emergent properties of the model (e.g. critical N content and self-thinning laws) and the needs for keeping this degree of model complexity have been tested. Further applications of the model are currently being developed including in silico experiments on the biodiversity-ecosystem functioning relationships and plant trait space analyses (Gross et al., in press). Further work will also be needed to extend the model to other plant functional types such as forbs and fully account for the effects of water deficit. This will open the possibility to model mechanistically multi-species grassland systems and help to understand the mechanistic links between individual, community and ecosystem responses to global changes.

Fig. 1 .

 1 Fig. 1. Unified Model Language (UML) diagram of Gemini model and its modules Canopt (vegetation) and Soilopt (soil, not described here).

Fig. 2 .

 2 Fig. 2. Algorithmic structure of the vegetation module (Canopt) in the Gemini model.

Fig. 3 .

 3 Fig. 3.(A) Simulated versus measured above-ground vegetative dry matter yields (DMY) of 13 grass species and of a Lolium cultivar grown in monocultures with four management treatments (C-, C+: 3 and 6 cuts per year; N-, N+: 120 and 360 kgN ha -1 yr -1 N fertilizer supply). Model was calibrated and optimized with data from the C-N+ treatment only. A restart procedure after spin-up runs was used to initialize simulations. (B) Simulated versus measured relative abundance of grass species grown in mixtures with four management treatments. Three mixtures of six species each were drawn at random from a pool of 12 species: (i) Dg, Fa, Fr, Lp, Pp, C. cristatus; (ii) Ap, Ao, Ae, Er, Hl, Tf; (iii) Ae, Dg, Er, Fa, Fr, Hl. The simulations started with a total initial mixture density of 1000 axes m -2 , with each species having the same initial axis density. See Table2for species abbreviations.

Table 1

 1

Table 2 .

 2

	Plant functional trait	Abbreviation	Unit	Linked reference	Model parameter	Model mechanism
	Leaf morphological trait
	Leaf length	LL	cm	[3,4,5]	L leaf0 , F sheath , b AL , a AL , bWL, aWL	Shoot morphogenesis
	Sheath length	SL	cm	[3]	F sheath	Shoot morphogenesis
	Leaf area	LA	cm 2	[1]	b AL , a AL	Shoot morphogenesis
	Leaf fresh mass	LFM	mg	[5]	bWL, aWL	Shoot morphogenesis
	Leaf dry mass	LDM	mg	[5]	bWL, aWL	Shoot morphogenesis
	Leaf dry matter content	LDMC	mgDM g -1 FM	[1,3,5]	LDMC	Shoot morphogenesis
	Specific leaf area	SLA	m 2 kg -1	[1-5]	bWL	Shoot morphogenesis
	Number of growing leaves	NG	nb tiller -1	[4]	ng leaf	Shoot morphogenesis
	Number of mature leaves	NM	nb tiller -1	[4]	nm leaf	Shoot morphogenesis
	Root morphological trait
	Primary root length	RL	cm	[5]	LL plast , Lroot1	Root morphogenesis
	Root maximal order	RO	nb	[5]	LL plast , rr plast , ordermax	Root morphogenesis
	Root lifespan	RLS	Degree day	[5]	Troot g1 , Troot m1	Root morphogenesis
	Primary root diameter	RD	mm	[1,5]	rr plast , rrroot1	Root morphogenesis
	Type of root reserve organ	TRO	Rhizome, stolon. . .	[1]	WrrCmax, WrrNmax	Root morphogenesis
	Plant morphological trait
	Vegetative elongated plant height	VE	cm	[1,5]	C	Shoot morphogenesis
	Tiller density	TD	tillers m -2	[4]	intcl, Tsen0	Population dynamic
	Plant phenological trait
	Earliness of growth	EG	-	[4]	bLER	Population dynamic
	Phyllochron	PH	Degree day	[2,4]	ph0	Shoot morphogenesis,
						Population dynamic
	Leaf and root physiological trait
	Leaf N resorption rate	RE	%	[2]	fns, RN, RP	Physiology
	Leaf mass loss at senescence	Mass loss	%	[2]	fcs, RC	Physiology
	Root N uptake capacity	Imax

  Reynolds,

	Environment	1	0..1	Canopt
	Manage ment	0..3	Gemini
			0..1	Soilo pt
	Animal	Cut	Fertilisation

TERMINATE SIMULATI ON PART 1 SIMU LATION START •

  Const ruct the model from sub-models: o 0 to n soil la yers, if 0 la yer an inert substrate is ass umed o 0 to n plant popul ations (grass or legu me) o Explicit or implicit morpholog y of shoots and roots , o Type of envi ronment: constant or va riable (dail y val ues from inpu t file ) Update de mograp hic statu s of lea ves an d roots: growin g, matu re or senescent . • Update organ number follo wing ne w emissions of leaves and roots • Compute the pote ntial daily extension of lea ves and roots, compute ass imilate de mand • Compute sene sced area of each organ and total senesced area. Compute fraction of senescence for flux calculations.

	o	Disturbance: None, Cutting, Grazing
	o	N supp ly type: dail y or afte r cuts
	• Impo rt dail y input data:
	o	Mea n temperatur e, temperatur e range, photoperiod, PPFD, precipitation, air moistur e	STAGE 2: COMPUTE SU BSTRATE SUPPL Y AND DEMAND
	o • Main options: Daily N input, cut or no cut Resta rt: save and/o r loa d 1 to n pop ulati ons o o Canopy la yer width	Supply • Photos ynthesis • Sene scence and re sorption • Soil N uptake and N fixation	D	i r • Maintenance o i r p y b d n a m e • Growth • Storage	t	y	o	d r	e	r	Defoliation event s • Export C and N substrates in proportion to defoliation height
				• Reserve remobili sation
				STAGE 3: PARTITION SU BST RATES MAT CHING SUPPLY AND DEMA ND
				•	In each canopy layer, compute co-limited leaf N contents per unit area
	PART 2 MA IN LOOP (Runge Kutta 4 or E ule r integration method )	•	Compute photos ynyhetic N de mand from the diff erence bet wee n act ual and co-limite d
					N lea f content
	1. Update resourc e bioavailability in each layer (soil and canopy) • Defoliation event s • Compute canop y layers number and in each layer leaf area and angle • Compute radiative balance in each layer	• •	Calculate the partitioning variable Q between leaf struc tures and leaf prot eins fro m the ratio between actual and co-limited lea f N contents In ca se of low substrate N availa bilit y de-repress root activity, reduce Q value an d hence partitioning to leaf proteins,
	•	Compute N balance in each soil layer	•	Compute functio nal balance betwe en shoo ts and roots and calculate P partitioning
					variable between shoot st ructu res and roots
	2. Update resourc e capture and partit ioning for each population • Morpho logical su b-model: o Tillering, leaf and root emission Biogeoche mical sub-model: •	• •	From P and Q, co mpute subst rate potential suppl y to the 3 st ructu ral compartments , Partition sub strate s to the 3 structural compartment s as the minimum of suppl y and of morphogenetical de mand .
		o	Compute photos ynthesis and N captu re
		o	Partition C and N substrates to maintenance , reserv es and grow th.
	•	Demography: Recruit new mat ure tillers
	3. Update object Population

• Morphological su b-model: dead and cut organs removed from axis, • If an y state va riable is belo w ze ro, population si mulati on stops and corresponding C and N is tran sferr ed to corr esponding soil la yer. • When no populati on remains, stop si mulation • Whe n simulation ti me is ov er, stop simulation . STAGE 1: UPDATE MOR PHO LOGY AN D DEMOGR APHY • STAGE 4: GROWTH • Compute actual growth for each growing lea f or root • Compute corresponding elo ngation of each gro wing leaf or roo t STAGE 5: TILLER RECRUITMENT • Update C and N subst rate and rese rve pools • Calculate the lowest canop y laye r which confe rs autot roph y to a daughte r axi s • Compute til ler recrui tment based on the ratio of su bstrate demand for autotrophy an d sup ply

Table 2

 2 Calibration of the vegetation module (Canopt) of Gemini. (A) Species specific parameter values for the 13 native grasses and for the control Lolium perenne (Clerpin) cultivar. (B) Initial value of the state variables. aAL, bAL, aWL, bWL leaf allometry parameters fitted from relationships between lamina area or volume and leaf length; bLER, power coefficient of temperature response of leaf elongation; C, parameter controlling the leaf lamina angle for a given canopy height and density; Fsheath, fraction of sheath in total leaf length; Lleaf0, maximal leaf length; LDMC, leaf lamina dry-matter per unit fresh-matter; ngleaf, nmleaf, growing and mature leaf number per tiller; ph0, phyllochron; LLplast, branching order dependency of potential root length; Lrootmax1, potential length of primary roots; NmRoots1, average number of primary roots; Orl, branching order dependency; RDMC, root dry-matter content; rrroot1, radius of primary roots; rrplast, root order dependency; Trootg1 and Trootm1, duration of primary growing and mature roots; Tr, thermal time interval between two successive root emission events; intcl, clonal integration; Tsen0, lifespan of an axis; RP, RC, RN, N protein and dry-matter resorption rates during leaf senescence; Su0, root N uptake potential rate; WrrCmax and WrrNmax, potential size of slow C and N reserve pools; WS, WP and WR, structural compartments of shoot structures; leaf proteins and roots, respectively; D, number of axes per population and per square meter; WC and WN, mass of C and N substrate compartments; WrC and WrN, mass of C and N labile reserve compartments; WrrC and WrrN, mass of C and N stable reserve compartments. Species abbreviations are: Alopecurus pratensis (Ap), Anthoxanthum odoratum (Ao), Arrhenatherum elatius (Ae), Dactylis glomerata (Dg), Elytrigia repens (Er), Festuca arundinacea (Fa), Festuca rubra (Fr), Holcus lanatus (Hl), Lolium perenne (Lp), Lolium perenne cultivar (Cl), Phleum pratense (Php), Poa pratensis (Pp), Poa trivialis (Pt), Trisetum flavescens (Tf).

	Symbol	Ap	Ao	Ae	Dg	Er	Fa	Fr	Hl	Lp	Cl	Php	Pp	Pt	Tf
	(A) Full model
	Shoot morphology at axis scale
	a AL (cm 2 cm -1 ) a WL (cm 2 cm -1 )	0.138 0.00204	0.216 0.00340	0.074 0.00105	0.160 0.00278	0.054 0.00075	0.314 0.00610	0.143 0.00359	0.398 0.00572	0.050 0.00095	0.042 0.000828	0.055 0.00069	0.459 0.00718	0.043 0.00063	0.022 0.00030
	b AL (dimensionless)	1.158	1.110	1.358	1.201	1.472	1.081	0.913	0.927	1.359	1.387	1.419	0.747	1.311	1.736
	b LER (dimensionless)	2.08	1.65	1.9	1.8	1.9	1.29	2.1	1.15	1.42	2.00	1.2	1.35	2.23	2.12
	b WL (dimensionless)	1.158	1.110	1.358	1.201	1.472	1.081	0.913	0.927	1.359	1.387	1.419	0.747	1.311	1.736
	C (dimensionless)	0.12	0.16	0.12	0.12	0.12		0.08	0.12	0.12		0.10	0.10	0.18	0.12
	F sheath (0-1)	0.313	0.270	0.309	0.261	0.285	0.204	0.225	0.353	0.253	0.300	0.337	0.232	0.360	0.320
	L leaf0 (cm) LDMC (gDM g -1 FM)	57.0 0.215	31.6 0.192	43.1 0.182	52.0 0.193	48.9 0.205	54.0 0.215	29.3 0.250	49.4 0.174	40.9 0.185	41.6 0.166	33.8 0.206	35.3 0.232	21.7 0.236	32.4 0.211
	ng leaf (leaves axis -1 ) nm leaf (leaves axis -1 ) ph 0 ( • C d)	1.17 2.71 163.5	1.12 2.41 200.0	1.18 2.44 162.1	1.79 2.53 142.0	1.21 2.53 157.0	1.01 2.00 220.0	1.03 2.70 239.3	1.14 2.38 130	1.16 2.42 154.4	1.22 2.28 212.0	1.36 2.73 120.0	1.11 2.23 220.0	1.18 3.16 115.0	1.23 2.53 200.0
	Root morphology at axis scale
	LL plast (dimensionless)	1.114	0.589	0.691	0.764	0.866	1.500	0.752	0.607	0.837	0.691	0.626	0.599	0.541	0.462
	L rootmax1 (cm) Nm Roots1 (roots axis -1 )	52 6	39 8	42 6	50 8	49 6	52 8	32 6	50 6	39 6	48 6	33 8	45 6	32 8	38 6
	OrL (dimensionless)	3.50	2.20	3.30	3.05	3. 80	1.4	2.90	2.10	4.50	3.00	2.30	2.20	4.27	2.70
	RDMC (gDM g -1 FM)	0.130	0.110	0.170	0.144	0.147	0.139	0.121	0.136	0.118	0.110	0.170	0.150	0.120	0.130
	rr root1 (cm)	0.0285	0.0161	0.0215	0.0221	0.0227	0.0310	0.0210	0.0160	0.0238	0.0170	0.0157	0.0171	0.0135	0.0147
	rr plast (dimensionless) T rootg1 ( • C d) T rootm1 ( • C d) Tr ( • C d)	1.114 300 500 80	0.589 500 200 28	0.691 300 500 40	0.764 300 500 45	0.866 300 500 49	1.500 300 500 52	0.752 300 500 65	0.607 300 500 45	0.837 300 600 150	0.691 300 500 50	0.626 300 500 40	0.599 300 500 70	0.541 300 500 20	0.462 300 500 32
	Morphology at population scale
	Intcl (0-1)	0.185	0.22	0.22	0.255	0.205	0.245	0.258	0.25	0.22	0.18	0.21	0.23	0.15	0.24
	Tsen 0 ( • C -1 d -1 )	0.0051	0.0051	0.007	0.0074	0.0052	0.0043	0.00425	0.00525	0.0073	0.0043	0.007	0.0035	0.007	0.0055
	Chemical composition
	fcs (gC g -1 DM)	0.425	0.433	0.429	0.429	0.419	0.432	0.441	0.424	0.439	0.435	0.432	0.434	0.419	0.432
	fcr (gC g -1 DM)	0.414	0.433	0.428	0.435	0.428	0.422	0.427	0.425	0.428	0.428	0.437	0.425	0.417	0.425
	fns (gN g -1 DM)	0.0173	0.0173	0.0238	0.0094	0.0106	0.0087	0.0137	0.0092	0.0182	0.0099	0.0123	0.0137	0.0224	0.0193
	fnr (gN g -1 DM)	0.0112	0.01198	0.0124	0.0107	0.0129	0.0127	0.0133	0.0130	0.0149	0.0124	0.0124	0.0109	0.0166	0.0101
	Physiology (@ T ref = 20 • C)
	R C (0-1)	0.362	0.522	0.279	0.408	0.372	0.258	0.302	0.498	0.251	0.385	0.469	0.269	0.368	0.319
	R N (0-1)	0.362	0.522	0.279	0.408	0.372	0.258	0.302	0.498	0.251	0.385	0.469	0.269	0.368	0.319
	R P (0-1) Su 0 (gN m -2 root d -1 ) W rrCmax (gC axis -1 ) W rrNmax (gN axis -1 )	0.421 0.1318 0.01 0.001	0.592 0.0618 0.01 0.001	0.253 0.2376 0.01 0.001	0.623 0.1282 0.01 0.001	0.585 0.1810 0.02 0.002	0.546 0.0750 0.04 0.004	0.490 0.0928 0.01 0.001	0.627 0.0710 0.01 0.001	0.330 0.0564 0.01 0.001	0.635 0.0960 0 0	0.590 0.0485 0.01 0.001	0.453 0.1404 0.01 0.001	0.377 0.0933 0 0	0.398 0.0493 0.01 0.001
	(B) Initial values of state variables
	D (axes m -2 )	2444	2510	3514	3028	3079	1314	8474	3652	5270	5472	4645	4554	4846	5358
	W C (gC axis -1 ) W N (gN axis -1 ) W P (gDM axis -1 ) W R (gDM axis -1 ) W rC (gC axis -1 ) W rN (gN axis -1 ) W rrC (gC axis -1 ) W rrN (gN axis -1 ) W S (gDM axis -1 )	0.029 0.017 0.004 0.106 0.025 0.017 0.022 0.006 0.060	0.025 0.020 0.004 0.086 0.022 0.020 0.021 0.006 0.049	0.029 0.014 0.005 0.074 0.026 0.014 0.022 0.005 0.031	0.029 0.017 0.004 0.106 0.024 0.017 0.022 0.006 0.064	0.019 0.012 0.003 0.077 0.015 0.012 0.023 0.005 0.047	0.042 0.022 0.007 0.114 0.078 0.036 0.077 0.054 0.152	0.009 0.005 0.001 0.035 0.007 0.004 0.011 0.001 0.031	0.017 0.011 0.003 0.061 0.014 0.011 0.016 0.005 0.059	0.014 0.012 0.002 0.053 0.012 0.012 0.015 0.005 0.028	0.023 0.007 0.002 0.078 0.021 0.007 0.000 0.000 0.026	0.014 0.009 0.002 0.049 0.012 0.009 0.015 0.004 0.031	0.016 0.007 0.002 0.053 0.013 0.007 0.016 0.004 0.046	0.009 0.006 0.002 0.029 0.007 0.006 0.012 0.003 0.022	0.013 0.008 0.002 0.040 0.011 0.008 0.014 0.004 0.028
	Parameter abbreviations:

Table 3

 3 Statistical indices used for Gemini model evaluation comparing simulated and measured values of (A) above-ground dry matter yield of 13 monocultures and (B) species abundance in three different mixtures of six species, each being cultivated in four management treatments (C-N+, C+N+, C-N-and C+N-).

	Four complemen-

J.-F.Soussana et al. / Ecological Modelling 231 (2012) [134][135][136][137][138][139][140][141][142][143][144][145]

Acknowledgements

This study was financially supported by the French ANR DISCOVER (ANR-05-BDIV-010-01) and ANR QDIV (ANR-05-BDIV-009-01) projects and by the FEDER, 'l'Europe s'engage en région Auvergne'. V. Maire was funded by a Ph-D grant of French Research Ministry. We thank C. Picon-Cochard for constructive comments on root parameterization, Hoglah Dasari and two anonymous referees who provided helpful comments on previous versions of the manuscript.


Julia Dominguez Escobar Dominguez Escobar

Arnaud Chastanet

Alvaro H A H Crevenna

Vincent V Fromion

Roland R Wedlich-Söldner

Rut R Carballido Y Lopez

Julia Domínguez-Escobar

Rut Carballido-López

Processive Movement of MreB-Associated Cell Wall Biosynthetic Complexes in Bacteria

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

lighted examples of proximal opposing movements (Fig. 4, E and F; fig. S15; and movies S13, A to D). Thus, PGEM and MreB filaments move in both directions around the cell, and we could not resolve any coordination along the cell length.

It thus appears that the coupled motions of the PGEM and MreB reflect the active process of cell wall synthesis: a circumferential motion of disconnected MreB-PGEM complexes moving around the cell in both directions, synthesizing discrete radial bands of PG oriented perpendicular to the cell length (fig. S16). This model is consistent with the arrangement of B. subtilis PG observed by atomic force microscopy (19).

MreB filaments are required for elongationspecific PG synthesis (2,20), suggesting that they are integral components of these translocating machines. These filaments may serve as coordinating scaffolds to link the PGEM to the enzymes that synthesize PG precursors (1-4, 10) (fig. S17). We cannot completely rule out the contribution of polymer dynamics to these motions, because there are no methods to inhibit MreB polymerization without disrupting existing filaments. However, we did not observe any directed motion in the absence of PG synthesis, even with high-precision measurements (fig. S13), suggesting that PG synthesis is the predominant process driving these motions. If these motions are driven by MreB polymerization or another process, this would require induction of equivalent rigor states during depletions of all PGEM components and the antibiotic inhibition of PG cross-linking, polymerization, and precursor synthesis.

Rather than a contiguous helical structure, these observations reveal the mobile, fragmented nature of MreB. Thus, although MreB is required for rod-shape maintenance, it cannot function as a cell-spanning structure, much less a coherent "cytoskeleton" in B. subtilis. It remains to be determined how the short-range activities of these independent biosynthetic complexes impart a long-range order to the cell wall.

The peptidoglycan cell wall and the actin-like MreB cytoskeleton are major determinants of cell shape in rod-shaped bacteria. The prevailing model postulates that helical, membrane-associated MreB filaments organize elongation-specific peptidoglycan-synthesizing complexes along sidewalls. We used total internal reflection fluorescence microscopy to visualize the dynamic relation between MreB isoforms and cell wall synthesis in live Bacillus subtilis cells. During exponential growth, MreB proteins did not form helical structures. Instead, together with other morphogenetic factors, they assembled into discrete patches that moved processively along peripheral tracks perpendicular to the cell axis. Patch motility was largely powered by cell wall synthesis, and MreB polymers restricted diffusion of patch components in the membrane and oriented patch motion.

T he peptidoglycan (PG) layer, or sacculus, which is composed of glycan strands crosslinked by peptide bridges, forms a loadbearing network that maintains bacterial cell shape. Synthesis and chemical composition of PG are well understood (1), but the structure of the sacculus and the mechanisms controlling its growth remain elusive. It is currently assumed that the actin-like MreB proteins form filamentous helical structures along the membrane, which direct cell wall growth by positioning multienzyme complexes that mediate sidewall elongation (2). These elongation complexes are thought to contain the essential transmembrane proteins MreC and MreD, RodA and RodZ, PG hydrolases, and penicillin-binding proteins (PBPs), the enzymes that catalyze PG elongation and cross-linking (3,4). Studies using fluorescently labeled vancomycin (Van-FL) have revealed helical incorporation patterns of PG precursors into the sidewall in B. subtilis (5) and Caulobacter crescentus (6), but could not resolve the dynamics of cell wall synthesis. We used total internal reflection fluorescence microscopy (TIRFM), a sensitive method for studying events at cell surfaces (7,8), to observe the dynamics of the MreB cytoskeleton and its relationship with cell wall growth.

Functional green fluorescent protein (GFP) fusions to the three MreB isoforms in B. subtilis (MreB, Mbl, and MreBH) expressed at wildtype levels (fig. S1) formed discrete patches in exponentially growing cells (Fig. 1A). Patches were restricted to the cell periphery (fig. S2A) and exhibited continuous movement along linear tracks roughly perpendicular to the long axis of the cell (Fig. 1B and movie S1). To reconcile our findings with the helical structures described for MreB proteins (2, 9), we simultaneously imaged cells by TIRFM and conventional epifluorescence. Owing to the increased depth of field, MreB patterns visualized by epifluorescence could be misinterpreted as "helical" (Fig. 1C). In addition, MreB localized to transverse bands as cells entered stationary phase (fig. S2B). We found similar patch-like localization and dynamics of MreB in the Gram-negative bacteria Escherichia coli and C. crescentus (fig. S2C), suggesting that these are widely conserved features in bacteria.

Detailed analysis of patch dynamics using kymographs (10) showed that patches moved bidirectionally across the cell at constant velocity (fig. S2D). We often observed reversal and crossover of patches (fig. S2D), but did not see abrupt patch appearances or disappearances at midtrajectory, indicating tight association of MreBs with the cell periphery. To determine patch velocities, we measured the slopes of kymograph traces. Patches formed by MreB, Mbl, and MreBH moved at significantly different speeds (Fig. 1D, fig. S3, and table S1; P < 0.001 for all pairs). Speeds were mildly dependent on growth temperature (fig. S4A) but were not affected by high concentrations of magnesium (fig. S4B) or untagged endogenous copies of the respective MreB isoforms (fig. S4C and table S1). Patch trajectories in maximum projections were oriented at angles close to 90°r elative to the long axis of the cell (Fig. 1E and table S1). Although tracks were evenly spaced along the length of the cell, with major distance peaks between 0.5 and 1 mm (Fig. 1F), high variability of autocorrelation between cells (fig. S5, A to C) and over time (fig. S5D) argued against a constrained periodic structure. However, we did find a correlation between numbers of tracks and cell length (Fig. 1G), suggesting an average distance between tracks of ~0.5 mm, as previously reported (11).

MreB, Mbl, and MreBH have been reported to colocalize and interact extensively with each other (12,13). To determine whether differences in patch speeds (Fig. 1D) reflected distinct isoform composition or simply different strain backgrounds, we performed two-color TIRFM on pairs of GFP-and mRFPruby-(monomeric RFP, red fluorescent protein) tagged proteins. All tested pairs displayed extensive overlap (Fig. 1H) with colocalization in more than 75% of kymograph traces (Fig. 1I and table S3) and moved at similar speeds (fig. S1, and movie S3). All other PBPs and LytE localized to patches that randomly moved along the cell surface (fig. S7A). In cells lacking pbpH, PBP2a patches displayed directional movement along circumferential tracks (Fig. 1, J and K, table S1, and movie S3), consistent with the proposed reciprocal redundancy of these two transpeptidases in sidewall synthesis (17). Furthermore, the integral membrane protein RodA, which has been linked to cell wall elongation and to PbpH and PBP2a (17,18), also formed circumferentially moving patches (Fig. 1, J and K, table S1, and movie S3). Finally, the orientation of Van-FL-labeled tracks was also centered around 90°C (Fig. 1K) instead of helical (5). Whereas MreBs and MreC patches exclusively displayed circumferential motion, MreD, PBP2a, PbpH, and RodA patches also frequently exhibited rapid diffusion along the membrane (fig. S7B). Circumferentially motile PbpH and RodA strongly colocalized with MreB and Mbl (Fig. 1, L and M, fig. S6, and table S1). Thus, MreBs, MreC/D, RodA, and PbpH/2a form part of cell wall elongation complexes that display circumferential processive motility in B. subtilis cells. What is the molecular basis for the processive movement? One possibility would be treadmilling of MreB, as previously suggested in B. subtilis (11) and C. crescentus (19). However, several dynamic behaviors of MreB patches such as fusion and fission (Fig. 2A and movie S4), reversal of direction (Fig. 2A), or splitting (Fig. 2B) argue against treadmilling-driven motion. To directly test for turnover of MreBs within patches, we performed FRAP (fluorescence recovery after photobleaching) experiments. When individual GFP-MreB or GFP-Mbl patches were partially bleached, their motility was unaffected and no fluorescence recovery occurred during their movement across the cell (Fig. 2, C andD). Similarly, when a whole cell was bleached, with the exception of a single GFP-Mbl patch (Fig. 2E), the signal from this patch did not decrease during its movement across the cell, again ruling out treadmilling-driven motion. Absence of treadmilling does not, however, preclude polymerization and slow global turnover of MreB (20), which is predicted from the correlation between cell length and patch numbers (Fig. 1F).

G-MreB Mbl-R G-MreBH R-MreB G-MreBH Mbl-R G-Mbl Mbl-R Mbl-G R-Mbl % traces 100 0 G-PbpH R-MreB R-Mbl RodA-G R-MreB RodA-G R-Mbl G-PbpH
Alternatively, a motive force for MreB patches, and thus for the elongation complexes, could be provided by PG synthesis itself. To test this possibility, we treated cells with vancomycin (100 mg/ml) to inhibit incorporation of PG precursors into the sacculus (21), or phosphomycin (700 mg/ml), which blocks the first cytosolic step in PG precursor synthesis (22). Addition of either antibiotic completely stopped movement of PbpH, MreB, and Mbl patches (Fig. 3, A andB). After removal of the drugs, motility gradually resumed (Fig. 3, C and D, and movie S5). Finally, upon disruption of the PG backbone with lysozyme, patch movement also ceased (Fig. 3E and movie S6), and many patches rapidly diffused in the cytosol (Fig. 3E and movie S6). If PG assembly directly drives motility of MreB patches, the concentration of PG precursors available should influence patch speed. Reducing the precursor pool with low concentrations (1 to 5 mg/ml) of phosphomycin did not affect the direction of MreB and Mbl patches (fig. S9A), but their velocities were significantly lower (Fig. 3F and table S1). FRAP experiments showed no turnover of MreB or Mbl patches in the presence of either high or low concentrations of phosphomycin (fig. S8), confirming the absence of MreB turnover in patches. To target PG synthesis genetically without inducing drastic changes in cell morphology, we analyzed MreB, Mbl, and RodA patches in mutants null for pbpH and pbpA (encodes PBP2a), which display a wild-type growth rate and morphology during vegetative growth (17,23). Both deletions significantly slowed down MreB and Mbl patches (Fig. 3G), suggesting a role for transpeptidases in setting the pace of new PG strand assembly. Indeed, the essential transpeptidase PBP2 of E. coli plays a direct role in glycan strand synthesis apart from its role in cross-linking (24). Again, track orientation was unaffected (fig. S9B). Mutants in pbpH and pbpA were slightly impaired in growth (fig. S9C), consistent with a link between patch motility and cell growth.

Our results indicate that cell wall synthesis itself provides the driving force for the processive motility of sidewall elongation complexes. This is also supported by recent analysis showing that PG glycosyl transferases are processive enzymes (25). What, then, is the biological function of the essential actin-like MreB? MreB positions cytosolic enzymes involved in PG precursor synthesis in C. crescentus (6,14). Thus, MreB patches might simply act as passive scaffolds for synthesis and export of PG precursors and/or regulate their delivery to sites of active PG synthesis (24). Our observation that RodA and PBPs localize not only to MreB-associated motile elongation complexes but also to rapidly diffusing patches (fig. S7A) suggested an alternative scenario, in which MreB polymers might actively restrict and/or control the mobility of cell wall elongation complexes. To test this hypothesis, we monitored patch motility in the absence of individual MreB isoforms. Patches moved much faster in the mreB mutant, whereas deletion of mbl or mreBH had no significant effect (Fig. 4A, fig. S10, and table S1). In addition, GFP-Mbl patches in a DmblDmreB background exhibited less uniform directionality (Fig. 4B), sometimes even following trajectories along the cell axis (Fig. 4C and movie S7). An even stronger effect was observed on RodA and PbpH: Processively motile patches were almost completely elimi-nated, and the few remaining covered very short distances, albeit faster (Fig. 4, D andE). Thus, MreB may function directly in restricting and/or organizing motility of cell wall elongation complexes, although an additional role in the recruitment of PG precursors cannot be excluded.

We propose a model for sidewall elongation in B. subtilis, where motile membrane-associated elongation complexes insert new PG along bands largely perpendicular to the long cell axis. Old PG strands are used as guiding scaffolds, and motility of the complexes is powered by PG polymerization. MreB-like polymers restrict the diffusion of the complexes within the membrane to achieve processive and correctly oriented movement. A similar scenario has been proposed for plant cells, where processive motility of cellulose synthase complexes is suggested to be constrained by cortical microtubules (26).

  S6 and tableS1), indicating that the three MreB isoforms coexist in motile patches.MreB and Mbl have been implicated in the spatial organization of lateral cell wall synthetic

Fig. 1 .

 1 Fig. 1. Motile cell wall elongation complexes. (A) GFP-MreB, GFP-Mbl, and GFP-MreBH (left to right) imaged by TIRFM. Data in (A) to (G) are from strains 3723, 2523, and 2566J (table S4). (B) Movement of a GFP-MreB patch (arrow) and trajectory in maximum projection. Time is in seconds. (C) Simultaneous visualization of GFP-Mbl by TIRFM and epifluorescence (EPI). (Left) Corresponding bright-field (BF) image. (D) Typical kymograph and patch speed boxplots of MreB isoforms. (E) Typical maximum projection and angle boxplots for MreB isoform trajectories. (F) Distribution of GFP-MreB trajectories in maximum projection, linescan (along dotted line) and intensity correlation function (ICF). (G) Number of patch trajectories increases with cell length (linear fit: R 2 = 0.61). (H) Colocalization of GFP-MreB and Mbl-RFP in maximum projections, linescan (dotted line) and kymographs. (I) Quantification of MreBs colocalization from kymograph traces. G: GFP; R: RFP; gray bars: colocalization; green and red bars: single GFP and RFP color traces. (J) Boxplots of patch speed of morphogenetic proteins. (K) Orientation of patch trajectories of morphogenetic proteins and Van-FL bands. (L) Colocalization of GFP-PbpH and RFP-MreB in maximum projections, linescan, and kymograph traces. (M) Quantification of MreB and Mbl colocalization with PbpH and RodA from kymograph traces. Labels as in (I). Scale bars: 1 mm.

Fig. 2 .

 2 Fig. 2. Patch motility is not driven by treadmilling. (A) Time series and kymograph showing fission and fusion (arrowhead) and reversal of a GFP-Mbl patch (asterisk, RWSB10). Red: TIRFM; green: epifluorescence; blue (cell outline): bright field. (B) Time series and kymograph (dotted line) showing splitting of a GFP-MreB patch (RWSB1). (C) TIRF-FRAP with partial photobleaching of a moving GFP-MreB patch (RWSB1). (D) Kymographs along patch traces of GFP-Mbl in the DmrebDmbl background (RWSB10), with corresponding intensity profiles showing lack of fluorescence recovery upon partial bleaching. (E) Inverse FRAP (iFRAP). A cell (RWSB10) was bleached within the region outlined, omitting only a GFP-Mbl patch. Kymograph (dotted line) shows movement of the patch with no loss of fluorescence. Scale bars: 1 mm; time bars: 30 s in (A) and (B); 10 s in (D) and (E). Time is in seconds. Asterisks and arrow at initial position in (C) and (E).A B

Fig. 3 .

 3 Fig. 3. Contribution of PG synthesis to patch motility. (A and B) Immobile MreB, Mbl, and PbpH patches after treatment with vancomycin (100 mg/ml) for 8 min (A) and with phosphomycin (700 mg/ml) for 30 min (B). Time is in seconds. (C and D) Effects of vancomycin (C) and phosphomycin (D) are reversible. Kymographs were taken before (-) and immediately after addition of the drug (0) and at the indicated times (in minutes) after washout. Red arrows: partial recovery of motility. (E) Kymographs of GFP-Mbl patches (2523) showing partial arrest (a), complete arrest (b), or diffusive motility (c) after treatment with lysozyme (1 mg/ml) for 5 min. (F and G) Patch speed boxplots upon low phosphomycin (1 to 5 mg/ml) treatment relative to untreated cells (F) and in the presence and absence of pbpH or pbpA (G). Scale bars: 1 mm; time bars: 30 s.

Fig. 4 .

 4 Fig. 4. Roles of MreBs in patch motility. (A and B) Boxplots of patch speeds (A) and trajectory angles (B) in the presence and absence of mreB. (C) Movement of a GFP-Mbl patch (asterisk) along the cell axis in the DmrebDmbl strain (RWSB10). Red: TIRFM; green: epifluorescence; blue (cell outline): bright field. Time is in seconds. (D) Normalized PbpH and RodA patch motility (patch trajectories per kymograph trace per minute) after deletion of mreB or mbl. (E) Typical kymographs for PbpH and RodA patches in wildtype and mreB backgrounds. Scale bars: 1 mm; time bar: 30 s.A

Max Planck Institute of Biochemistry, Am Klopferspitz 18, D-82152 Martinsried, Germany.

INRA, UMR1319 Micalis, F-78352 Jouy-en-Josas, France.

AgroParisTech, UMR Micalis, F-78350 Jouy-en-Josas, France.

INRA, UR1077 Mathématique, Informatique et Génomes, F-78352 Jouy-en-Josas, France. *These authors contributed equally to this work. †To whom correspondence should be addressed. E-mail: wedlich@biochem.mpg.de (R.W.-S.); rut.carballido-lopez@jouy. inra.fr (R.C.-L.) www.sciencemag.org SCIENCE VOL 333 8 JULY 2011

on January 4, 2013 www.sciencemag.org Downloaded from

Acknowledgments: We thank P. Graumman, J. Errington, R. Carballido-Lopez, D. Kearns, D. Popham, and D. Scheffers for strains; D. Kahne, S. Walker, T. Bernhardt, and T. Uehara for discussions; Q. Justman for editing; S. Layer for inspiration; T. Emonet for MicrobeTracker; and the HMS Nikon Imaging Center for the use of their facilites. Funding was provided by NIH grants R01-GM073831 (D.Z.R.), R01-GM096450 (X.Z.), and R01-GM39565 (T.M.). X.Z. is a Howard Hughes Medical Investigator. E.C.G. was supported by the Helen Hay Whitney Foundation. Acknowledgments: We thank R. Daniel and L. J. Wu for sharing unpublished strains and for helpful suggestions, D.-J. Scheffers for strains, N. Campo for the antibody against GFP, and P. Hardy for manuscript editing. This work was supported by a Human Frontier Science Program (HFSP) Young Investigator grant (R.C.-L. and R.W.-S., HFSP-RGY0067/2009-C); a European Molecular Biology Organization short-term fellowship (R.C.-L), The French National Research Agency (ANR) (R.C.-L., ANR-08-JCJC-0024-01), and the Max Planck society (R.W.-S.).

Vladimir V. Rogov, 4,5 Nathan R. Brady, 6 Benjamin Richter, 1 Jelena Korac, 7 Oliver Waidmann, 1 Chunaram Choudhary, 3 Volker Dötsch, 4 Dirk Bumann, 2 Ivan Dikic 1,7 * Selective autophagy can be mediated via receptor molecules that link specific cargoes to the autophagosomal membranes decorated by ubiquitin-like microtubule-associated protein light chain 3 (LC3) modifiers. Although several autophagy receptors have been identified, little is known about mechanisms controlling their functions in vivo. In this work, we found that phosphorylation of an autophagy receptor, optineurin, promoted selective autophagy of ubiquitin-coated cytosolic Salmonella enterica. The protein kinase TANK binding kinase 1 (TBK1) phosphorylated optineurin on serine-177, enhancing LC3 binding affinity and autophagic clearance of cytosolic Salmonella. Conversely, ubiquitin-or LC3-binding optineurin mutants and silencing of optineurin or TBK1 impaired Salmonella autophagy, resulting in increased intracellular bacterial proliferation. We propose that phosphorylation of autophagy receptors might be a general mechanism for regulation of cargo-selective autophagy.

M acroautophagy (hereafter referred to as autophagy) is an evolutionarily conserved catabolic process by which cells deliver bulk cytosolic components for degradation to the lysosome (1-4). Selectivity in cargo targeting is mediated via autophagy receptors that simul-taneously bind cargoes and autophagy modifiers, autophagy-related protein 8 (ATG8)/ microtubuleassociated protein light chain 3 (LC3)/g-aminobutyric acid receptor-associated protein (GABARAP) proteins, which are conjugated to the autophagosomal membranes (5,6). The regulatory mechanisms controlling the spatiotemporal dynamics of the autophagy receptor-target interaction in cells remain unclear (7).

Frankfurt Institute for Molecular Life Sciences and Institute of Biochemistry II, Goethe University School of Medicine, Theodor-Stern-Kai 7, D-60590 Frankfurt (Main), Germany. 2 Infection Biology, Biozentrum, University Basel, Klingelbergstr. 50/70, CH-4056 Basel, Switzerland. 3 The Novo Nordisk Foundation Center for Protein Research, Faculty of Health Sciences, University of Copenhagen, Blegdamsvej 3, 2200 Copenhagen, Denmark. 4 Institute of Biophysical Chemistry and Center for Biomolecular Magnetic Resonance, Goethe University, Marie Curie Strasse 9, D-60439 Frankfurt (Main), Germany. 5 Institute of Protein Research, 142290, Pushchino, Russia. 6