Thomas Koffel
email: koffelth@msu.edu.

Simon Boudsocq

Nicolas Loeuille

Tanguy

Tanguy Daufresne

Facilitation-versus competition-driven

Keywords: Nitrogen fixation, facilitation, zero net growth isocline, alternative stable states, primary succession, ecosystem development, resource-ratio theory

Introduction

Primary succession, the sequential replacement of plant species through time following the appearance of a bare substrate, is a central concept in ecology [START_REF] Walker | Primary Succession and Ecosystem Rehabilitation[END_REF][START_REF] Shugart | Ecological succession and community dynamics[END_REF], and occurs on virtually any substrate, from sand dunes [START_REF] Cowles | The physiographic ecology of Chicago and vicinity: a study of the origin, development, and classification of plant societies[END_REF][START_REF] Olson | Rates of succession and soil changes on southern Lake Michigan sand dunes[END_REF] to glacial moraines [START_REF] Chapin | Mechanisms of primary succession following deglaciation at Glacier Bay, Alaska[END_REF]. [START_REF] Odum | The strategy of ecosystem development[END_REF] defines succession through three characteristics: (i) it is a reasonably directional and orderly process, (ii) it results from the reciprocal interaction between the plant community and its physical environment, (iii) it culminates in a climax maximizing biomass and stability.

Nitrogen (N)-fixing organisms such as legumes and actinorhizal plants usually thrive during primary succession due to an initial lack of available N [START_REF] Walker | Primary Succession and Ecosystem Rehabilitation[END_REF]. In Glacier Bay, Alaska, succession usually includes a broad range of N-fixing species [START_REF] Walker | Nitrogen fixers and species replacements in primary succession[END_REF], from pioneers like 'black crust'-forming cyanobacteria [START_REF] Worley | The "black crust" phenomenon in Upper Glacier Bay, Alaska[END_REF][START_REF] Schmidt | The earliest stages of ecosystem succession in high-elevation (5000 metres above sea level), recently deglaciated soils[END_REF]) and lichen of the genus Stereocaulon, to mid-and late colonists like Drummond's aven (Dryas drummondii ; [START_REF] Lawrence | The role of Dryas drummondii in vegetation development following ice recession at Glacier Bay, Alaska, with special reference to its nitrogen fixation by root nodules[END_REF]) and green alder (Alnus viridis). Conversely, because of biomass turnover and recycling, fixed N accumulates in the soil and becomes available to the whole community [START_REF] Crocker | Soil development in relation to vegetation surface age at Glacier Bay, Alaska[END_REF][START_REF] Chapin | Mechanisms of primary succession following deglaciation at Glacier Bay, Alaska[END_REF][START_REF] Halvorson | Decomposition of lupine biomass by soil-microorganisms in developing Mount St-Helens pyroclastic soils[END_REF][START_REF] Crews | Organic matter and nitrogen accumulation and nitrogen fixation during early ecosystem development in Hawaii[END_REF]Walker et al. 2003;[START_REF] Kohls | An assessment of soil enrichment by actinorhizal N2 fixation using d15N values in a chronosequence of deglaciation at Glacier Bay, Alaska[END_REF]. Such facilitation mechanisms between plants [START_REF] Callaway | Positive interactions among plants[END_REF] are potential drivers of primary succession, as exemplified by the facilitation model of succession, or facilitation-driven succession, in which 'early colonizers modify the environment, thereby making it more suitable for later colonizers to invade and prosper and less suitable for early colonizers to survive' [START_REF] Connell | Mechanisms of succession in natural communities and their role in community stability and organization[END_REF].

In fact, N-fixing organisms also compete for the access to other essential resources such as light or phosphorus (P). Contrary to N availability, P availability in the pre-succession substrate is usually high, then decreases over time as it gets weathered from the bedrock and immobilized in the biotic compartments [START_REF] Walker | The fate of phosphorus during pedogenesis[END_REF][START_REF] Vitousek | Nutrient limitations to plant growth during primary succession in Hawaii Volcanoes National Park[END_REF][START_REF] Vitousek | Terrestrial phosphorus limitation: mechanisms, implications, and nitrogen-phosphorus interactions[END_REF]. These central yet contrasted roles played by N and P during succession have put them at the center of models of ecosystem development over the last decades [START_REF] Vitousek | Nutrient limitation and soil development: experimental test of a biogeochemical theory[END_REF][START_REF] Marleau | A stoichiometric model of early plant primary succession[END_REF][START_REF] Laliberté | Experimental assessment of nutrient limitation along a 2-million-year dune chronosequence in the south-western Australia biodiversity hotspot[END_REF][START_REF] Menge | Nitrogen and phosphorus limitation over long-term ecosystem development in terrestrial ecosystems[END_REF]. Both the absolute and relative richnesses of the abiotic environment in these two nutrients are expected to control succession, as formalized in the 'resource-ratio theory of succession' [START_REF] Tilman | The resource-ratio hypothesis of plant succession[END_REF]). Yet, it is still unclear how competition for P and facilitation for N interact during succession [START_REF] Callaway | Competition and facilitation: a synthetic approach to interactions in plant communities[END_REF], and when this results in facilitation-driven successions.

Models of primary succession have long been limited to verbal and conceptual approaches [START_REF] Clements | Plant Succession. An Analysis of the Development of Vegetation[END_REF][START_REF] Odum | The strategy of ecosystem development[END_REF][START_REF] Connell | Mechanisms of succession in natural communities and their role in community stability and organization[END_REF]. Their formalization first led to Markovian models of succession, i.e. based on phenomenological transitions between vegetated states, without specifying any competition mechanisms or making explicit the soil properties that mediates them [START_REF] Van Hulst | On the dynamics of vegetation: Markov chains as models of succession[END_REF][START_REF] Usher | Modelling ecological succession, with particular reference to Markovian models[END_REF]Siles et al. 2008). By contrast, population- [START_REF] Tilman | The resource-ratio hypothesis of plant succession[END_REF] or individualbased [START_REF] Huston | Plant succession: life history and competition[END_REF] resource competition models rely on the feedback between plants and their environment.

In the context of P acquisition and N fixation, recent resource competition models mostly consider short-scale transient dynamics including only a few different species [START_REF] Marleau | A stoichiometric model of early plant primary succession[END_REF] or a single undifferentiated plant compartment [START_REF] Menge | Nitrogen and phosphorus limitation over long-term ecosystem development in terrestrial ecosystems[END_REF]. Other succession models rely on a large number of species in a game theoretical perspective, but ignore N accumulation [START_REF] Tilman | The resource-ratio hypothesis of plant succession[END_REF][START_REF] Loreau | Ecosystem development explained by competition within and between material cycles[END_REF][START_REF] Sheffer | Biome-scale nitrogen fixation strategies selected by climatic constraints on nitrogen cycle[END_REF] or have focused on the succession end point rather than the whole developmental trajectory [START_REF] Menge | Evolutionary tradeoffs can select against nitrogen fixation and thereby maintain nitrogen limitation[END_REF][START_REF] Kylafis | Ecological and evolutionary consequences of niche construction for its agent[END_REF].

In this paper, we investigate autogenic succession by a community of N-fixers, differing in their abilities to fix N and acquire P. First, we use a mechanistic plant-soil model to describe how these species perform in the absence of other competing species, and classify them into three categories accordingly. Second, we embed this ecological model into a community assembly dynamics model to answer the following questions: how do competition and facilitation interact, and under what conditions does it lead to facilitation-driven succession? How does ecosystem stability change during succession? We conclude by summarizing these results into an enriched version of [START_REF] Tilman | The resource-ratio hypothesis of plant succession[END_REF] 'resource-ratio theory of succession'.

Ecological model and analysis

Plant-soil model

We consider a classical plant-soil model from resource competition theory [START_REF] Macarthur | Species packing and competitive equilibrium for many species[END_REF][START_REF] León | Competition between two species for two complementary or substitutable resources[END_REF][START_REF] Tilman | Resource Competition and Community Structure[END_REF][START_REF] Grover | Competition, herbivory and enrichment: nutrient-based models for edible and inedible plants[END_REF][START_REF] Chase | Ecological Niches: Linking Classical and Contemporary Approaches[END_REF] adapted to account for N-fixation, where a single species with population biomass B F grows on soil nitrogen (N) and phosphorus (P), with concentrations N and P respectively. N and P are supplied externally at constant influxes I N and I P and leave the system at leaching rates l N and l P . Following Leibig's law of the minimum, plant growth is limited by the least available nutrient:

g tot (P, N ) = min[g P (P ), g N (N ) + F ] (1)
where plant growth rate when P is limiting is given by g P (P ), and plant growth rate when N is limiting is the sum of growth on soil N denoted g N (N ) and maximal N-fixation rate F . Assuming constant plant stoichiometric coefficient q P and q N , plant growth translates into proportional P uptake equal to q P .g tot [START_REF] Tilman | Resource Competition and Community Structure[END_REF]. On the N side, we make the key assumption of facultative N-fixation, which is more realistic than purely obligate N-fixation (e.g. Tyrrell 1999), with important ecological consequences [START_REF] Agawin | Competition and facilitation between unicellular nitrogen-fixing cyanobacteria and non-nitrogen-fixing phytoplankton species[END_REF][START_REF] Menge | Facultative versus obligate nitrogen fixation strategies and their ecosystem consequences[END_REF]). More precisely, facultative N-fixation translates into soil N uptake being systematically privileged over atmospheric N-fixation, the latter only happening when soil N is limiting. This mechanism is justified by the substantial costs associated with N-fixation [START_REF] Gutschick | Evolved strategies in nitrogen acquisition by plants[END_REF]) and implicitly assumes that the plants can freely down-regulate their N-fixing activities. As a result, soil N uptake is equal to q N . min[g P (P ), g N (N )]. The interplay between P limitation and N-fixation leads to three different uptake regimes for the plant: P-limitation without N-fixation, P-limitation with partial N-fixation and N-limitation with maximal N-fixation (Fig 1 andAppendices A, B). Finally, plant biomass is lost at constant per capita mortality rate m and partially recycled to P and N with efficiency λ. The plant-soil population dynamics is then given by the following ordinary differential

dB F dt = [min[g P (P ), g N (N ) + F ] -m] B F (2a) dP dt = I P -l P P -q P min[g P (P ), g N (N ) + F ]B F + q P λmB F (2b) dN dt = I N -l N N -q N min[g P (P ), g N (N )]B F + q N λmB F (2c)
These equations are a simplified version of the model of [START_REF] Menge | Facultative versus obligate nitrogen fixation strategies and their ecosystem consequences[END_REF], not explicitly including organic nutrient pools, as they have no influence on our results at the ecological equilibrium. In practice, the associated extra parameters can be thought of being lumped in λ.

Regional species pool

We consider that the various N-fixing plant species that can establish during succession form a regional pool of strategies with different maximal fixation rates F , from plants with low fixing abilities that mostly rely on soil nitrogen to plants that can be completely independent from soil N referred to as "total fixers" [START_REF] Holter | N2(C2H2)-fixation in early stages of a primary succession on a reclaimed salt marsh[END_REF]. We make the key assumption that plants' maximal fixation rates F trade off against their ability to acquire P (see Table 1 for details, Appendix D for loosening this assumption, and Discussion). This is justified by the substantial costs associated with plants' specific adaptations to low-P soil conditions [START_REF] Zhu | The contribution of lateral rooting to phosphorus acquisition efficiency in maize (Zea mays) seedlings[END_REF][START_REF] Lynch | Rhizoeconomics: Carbon costs of phosphorus acquisition[END_REF], such as root morphology and architecture (e.g. lateral and cluster roots), root exudates, rhizosphere acidification and mycorrhizal symbioses, which leads to trade-offs between P-acquisition and other costly functions such as N-fixing ability, e.g. root-nodule allocation trade-offs [START_REF] Thuynsma | Phosphorus deficiency affects the allocation of below-ground resources to combined cluster roots and nodules in Lupinus albus[END_REF].

Graphical analysis and N-fixers classification

Before addressing succession dynamics, we first visualize the equilibrium properties of an ecosystem dominated by a given plant species with maximal fixation rate F along N and P external supplies. Doing so, we classify N-fixing species from the species pool into three categories -low, high and total N-fixers-according to their qualitative behavior, highlighting the different roles the N-fixing species will play during succession. Our model is adapted from the essential resource case of competition theory [START_REF] Tilman | Resource Competition and Community Structure[END_REF][START_REF] Grover | Resource Competition[END_REF][START_REF] Chase | Ecological Niches: Linking Classical and Contemporary Approaches[END_REF] with the inclusion of nitrogen fixing and nutrient recycling [START_REF] Daufresne | Plant coexistence depends on ecosystem nutrient cycles: extension of the resourceratio theory[END_REF]. It can therefore be studied graphically with the help of three ingredients: supply points, Zero Net Growth Isoclines (ZNGI) and impact vectors.

The supply point (S P , S N ) corresponds to N and P levels at equilibrium in the absence of any plant. It characterizes the baseline fertility of the bare substrate. Setting eq. (2b-2c) to zero with B = 0 gives S P = I P /l P and S N = I N /l N .

The supply point can be displayed in the resource plane (Black dot, Fig. 2b-d).

The ZNGI of a N-fixing species with biomass B F can be obtained by setting the growth eq. (2a) to zero. It gives the nutrient levels required for the plant to persist at equilibrium. The supply points situated above the ZNGI are the bare substrate conditions under which this species can establish itself starting from a very low density. As in the standard essential resource case [START_REF] Tilman | Resource Competition and Community Structure[END_REF], the ZNGI is a right-angle line whose corner's coordinates are Tilman's R * for N and P, the minimal resource required for growth (Fig. 2). We denote them P * and N * f , the latter subscript f distinguishing N * f from the plant's minimal requirements in the absence of fixation N * nf (See Appendix B for details). As N * nf > N * f , N-fixation enables the plant to grow under harsher N-limiting conditions, graphically lowering the position of the ZNGI corner. Because of the trade-off between maximal fixation rate F and P acquisition, P * and N * f are negatively correlated.

This can be visualized graphically when a large number of species with varying maximal fixation rates are sampled from the species pool and their ZNGIs displayed (Fig. 2a): all the ZNGIs' corner form a decreasing curve called the envelope [START_REF] Klausmeier | A model of flexible uptake of two essential resources[END_REF][START_REF] Danger | Does Liebig's law of the minimum scale up from species to communities?[END_REF][START_REF] Koffel | Geometrical envelopes: extending graphical contemporary niche theory to communities and eco-evolutionary dynamics[END_REF]. This continuum of strategies can be divided into three categories: low, high and total N-fixers (respectively Fig. 2b-d). Total N-fixers correspond to strategies with F > m, as they have N * f = 0. They can invade the bare substrate even in the complete absence of N supply, but are the poorest P

competitors (Fig 2c).

To understand the difference between low and high fixers, we need to introduce the impact vectors. The impact vectors describe the net impact of the plant population on the two nutrients, through uptake and recycling. Their coordinates can be read as the factors in front of B F in eq. (2b,2c). They are displayed on Fig. 2 for a few equilibrium points located on the ZNGI. Contrary to the low N-fixer case (Fig. 2b), the high and total N-fixer cases (Fig. 2c,d) present a soil N threshold denoted N * c under which plant presence has a net positive impact on N, with biotic supply by fixation and recycling exceeding soil uptake. This facilitation by and between N-fixing individuals under low N conditions leads to N accumulation and destabilizes the N-limited equilibria (Dashed ZNGI lower branch in Fig. 2c, compare with Fig. 2b).

In the high N-fixer case, facilitation leads to alternative stable states, namely, the non-zero plant population state and the bare substrate (Purple region in Fig. 2c; mathematical details in Appendix B). Indeed, N supply under the ZNGI is too low for N-fixers to invade from low density, but a larger population can persist under these conditions thanks to facilitation. This alternative stable state regime disappears for total N-fixers, as they can invade for any N availability (Fig. 2d). The presence of alternative stable states means that the vegetated ecosystem can undergo a sudden transition to the bare substrate [START_REF] Scheffer | Catastrophic shifts in ecosystems[END_REF][START_REF] Kéfi | Evolution of local facilitation in arid ecosystems[END_REF]. This sudden transition could be triggered by a slight change in environmental conditions such as a decrease in external phosphorous supply, or happen after a perturbation, e.g. fire or herbivore outbreak, that would be large enough to drive the N-fixing population under its recovery threshold.

Reading Fig. 2 as bifurcation diagrams along the nutrient supplies summarizes this analysis, with the white region corresponding to the bare substrate only (the plant can neither invade nor persist), the blue region the plant population only (the plant can invade and persist) and the purple region alternative stable states between the two (the plant cannot invade, but can persist).

Community assembly dynamics

Given the species pool described above, we now investigate the succession dynamics occurring among these species. Indeed, we suppose that these N-fixing species form a regional species pool [START_REF] Drake | The mechanics of community assembly and succession[END_REF][START_REF] Morton | Regional species pools and the assembly of local ecological communities[END_REF][START_REF] Chase | Ecological Niches: Linking Classical and Contemporary Approaches[END_REF], from which colonization attempts happen at random with the introduction of a few invading individuals into the succession site. To make the problem tractable, we assume a time scale separation between rare colonization attempts from the N-fixing species pool and faster ecosystem relaxation between these invasions. As we are interested in the role played by facilitation during succession, we also assume for the next two subsections that N availability is low enough for N-fixation to induce a positive feedback loop (S N < N * c ).

Successional step

Let us describe what happens when a focal invading species is introduced in the system formerly at equilibrium. This pre-invasion equilibrium situation corresponds either to the initial bare substrate, or follows the stabilization of a resident vegetated ecosystem after a previous successful colonization, as described in the previous section. In both cases, the preinvasion soil properties are given by their equilibrium nutrient concentrations, noted N and P , and completely determine the success of the invader. Denoting F and B F the maximal fixation rate and biomass of this invader, its success is quantified by the sign of its per capita growth rate when rare, the invasion fitness ρ [START_REF] Metz | How should we define 'fitness' for general ecological scenarios?[END_REF][START_REF] Geritz | Dynamics of adaptation and evolutionary branching[END_REF][START_REF] Geritz | Evolutionarily singular strategies and the adaptive growth and branching of the evolutionary tree[END_REF][START_REF] Hui | Defining invasiveness and invasibility in ecological networks[END_REF]):

ρ(F , P , N ) = 1 B F dB F dt P , N ;B F →0 + = min[g P ( P ), g N ( N ) + F ] -m (3)
where g P denotes phosphorus uptake by the invader and trades off against F . According to the previous section, if the pre-invasion ecosystem is the bare substrate, N and P are simply equal to the supplies S N and S P . Otherwise, N and P are fully determined by the equilibrium of the pre-invasion plant-soil model, i.e. located on the ZNGI of the resident plant species, whose maximal fixing rate is denoted F (like in Fig. 2b). If ρ(F , P , N ) is negative, invasion fails as the invader goes extinct and the resident ecosystem stays the same. If ρ(F , P , N ) is positive, invasion is successful and the invader's population gets large enough to impact the resident equilibrium. From there, several outcomes from this competition phase are technically possible when invading a resident population: 1) the invader becomes the new resident after excluding the former resident; 2) the invader and the resident reach an equilibrium where they coexist; or 3) invasion by the invader leads to a collapse of the whole ecosystem towards the bare substrate state, a particular case of a phenomenon known as 'the resident strikes back' [START_REF] Mylius | The resident strikes back: invader-induced switching of resident attractor[END_REF][START_REF] Geritz | Invasion dynamics and attractor inheritance[END_REF]. Under the assumption of low N availability (S N < N * c ), the positive feedback loop ensures that all strategies can only be P-limited at equilibrium, which, in virtue of Tilman's (1982) R * rule, excludes option 2, i.e. coexistence. As we also checked numerically that option 3, i.e. invasiondriven ecosystem collapse, does not happen, the only outcome here is option 1, i.e. exclusion of the resident by the invader.

This means that the community is always dominated by a single species throughout the succession.

As the invader replaces and becomes the new resident, it sets the nutrient levels at equilibrium on its own ZNGI.

This new equilibrium can have the bare substrate as an alternative stable state if strategy F cannot invade on its own (e.g. situation with high N-fixer depicted in Fig. 2c). When this happens, the presence of the resident F facilitated the establishment of F , by making the invasion of the latter possible [START_REF] Gerla | Photoinhibition and the assembly of light-limited phytoplankton communities[END_REF][START_REF] Rapaport | Some non-intuitive properties of simple extensions of the chemostat model[END_REF]. We call such an invader-resident replacement a facilitation-driven successional step. This facilitation-driven successional step is illustrated in Fig. 3 with two different strategies from the regional species pool, F 1 close to total N-fixation but not very competitive for P and F 2 with high N-fixation but more balanced N and P competitive abilities (F 1 > F 2 ). The supply point is such that the initial bare substrate is rich in P but highly N depleted, meaning that only F 1 , the pioneer species, can at first invade (S in Fig. 3a). However, the establishment of F 1 increases N availability at equilibrium due to positive impact (R 1 in Fig. 3). This leads to N level that are above F 2 requirements, making invasion by F 2 possible. Being better at P-acquisition, F 2 competitively excludes F 1 by pulling P levels down as it invades, and further increases N levels in the system (R 2 in Fig. 3). As a new resident, F 2 can persist after outcompeting F 1 thanks to the positive feedback loop on N, as S is located in F 2 alternative stable state region.

This successional step can be visualized using a graphical tool adapted from the eco-evolutionary literature, the Pairwise Invasibility Plot (PIP) [START_REF] Geritz | Dynamics of adaptation and evolutionary branching[END_REF]. It consists in plotting the sign of the invasion growth rate ρ against both the invader trait F on the y-axis, and the pre-invasion state (resident F or bare substrate) on the x-axis (Fig. 3b). The figure obtained is divided in two parts: the PIP proper on the left displaying invasion success of an invader F against a resident F , and an additional bar plot on the right displaying the invader success in the bare substrate (ρ > 0 in blue, ρ < 0 in white and gray Fig. 3b). The strategies with positive invasion on the bar plot correspond to the pioneer species.

Information on the status of the resident state is displayed on the plot, indicating whether this strategy cannot be resident for this supply point (black), or can be resident with (grey) or without an alternative stable state (white). The trajectory on the PIP illustrates the substitution sequence leading to the successional step.

Facilitation-driven succession and ecosystem development

The facilitation-driven successional step described in the previous section can be iterated by repeated successful colonizations events, leading to facilitation-driven successions. We assume equal probability of introduction for all the N-fixing species, i.e. no colonization-competition trade-off [START_REF] Tilman | Constraints and tradeoffs: toward a predictive theory of competition and succession[END_REF]). As discussed later, we expect facilitation-driven succession to be robust to the release of this assumption. Some general properties of this succession can be read directly on the PIP or deduced from the equilibrium properties of the plant-soil model. We also randomly generate successional trajectories and track the three quartiles of the distribution of probability to be in a given succession state through time (Fig. 4c-f), using the information contained in the PIP to assess the invasion success of invaders repeatedly chosen at random from the species pool.

In terms of community assembly, succession always starts with a pioneer species, i.e. a plant presenting a high enough maximal fixation rate to invade the bare substrate. This species can be displaced by another pioneer species, or facilitate its replacement by species with higher P-acquisition and lower fixation rates unable to invade by their own, themselves facilitating species with even lower fixation rates, etc, in an orderly monotonic process. However, as fixation rates decrease along succession, a decreasing range of strategies are able to invade, up to a certain strategy where no further invasion is possible (Fig. 4b). Such successional endpoint is located at the edge of the positive invasion region (blue, Fig. 4b), and is called an Evolutionary Stable Strategy (ESS) in the eco-evolutionary literature. Interestingly, this ESS is never attained but only approached asymptotically, as it is unable to invade any other strategy [START_REF] Geritz | Evolutionarily singular strategies and the adaptive growth and branching of the evolutionary tree[END_REF]. This leads to a progressive slowdown in late succession. As the ESS is located on the edge of the bistable region (gray, Fig. 4b), the late succession ecosystem is not globally stable, and a strong enough perturbation affecting plant biomass, such as fire or a sudden herbivore outbreak, could lead to the collapse of succession towards a non-vegetated state.

Such a succession of above-ground N-fixing strategies is associated with simultaneous below-ground changes in soil properties, leading to ecosystem development. Soil P starts at the initial substrate availability S P and then decreases, tending towards a reduced availability as we get close to the ESS (Fig. 4e). Conversely, soil N starts at S N and increases, asymptotically reaching a plateau (Fig. 4f). Such accumulation of soil N over time during succession seems at first contradictory as the strategies that establish in the course of succession are less and less efficient at N-fixation (Fig. 4c).

However, plant biomass also increases during succession, as plants better at exploiting P establish in the ecosystem (Fig. 4d). As total fixation is proportional to both N-fixer biomass and fixation rate, the increase in biomass compensates the decrease in per capita fixation rates, leading to an overall increase in fixation, and thus soil N, during succession. A last characteristic of facilitation-driven succession relies in the relative similarity between multiple trajectories despite the intrinsic stochastic nature of the colonization events (Fig. 4b-d). This is captured by the particular shape of the PIP, as the positive invasion region is restricted to a relatively narrow band.

A resource-ratio theory of succession

We now want to delimit exactly the nutrient supplies that lead to facilitation-driven succession, and explore other possible succession scenarios [START_REF] Connell | Mechanisms of succession in natural communities and their role in community stability and organization[END_REF]. Drawing such a bifurcation diagram along the nutrient supply gradients can be done easily using the 'envelope approach', a graphical tool adapting resource competition theory to community assembly dynamics [START_REF] Koffel | Geometrical envelopes: extending graphical contemporary niche theory to communities and eco-evolutionary dynamics[END_REF].

The results are summarized in Fig. 5. Four regions can be identified and labeled with corresponding PIPs. In the bare substrate region (0), N and P are so low that no strategy from the regional species pool can establish or persist. The facilitation-driven succession region (I) is found for low enough N (S N < N * c ) and high enough P supplies (on the right of the envelope). Interestingly, a small region (II) is located between (0) and (I). There, succession cannot start because no strategy can invade the bare substrate, but the community would go through succession and persist there if a population of high N-fixers were already present. This can be thought of as community-level alternative stable states between the bare substrate and a community that would undergo succession. Finally, region (III) emerges at high N supplies (S N < N * c ),

where fixation is too low for facilitation to happen. There, succession is purely competition-driven, as can be seen on the corresponding PIPs: most strategies can invade the initial bare substrate, including the ESS, and the establishment of any strategy only decreases the chances of further invasion by other strategies. Contrary to the facilitation-driven succession case, trajectories are not monotonous as the trait F can oscillate around the ESS, and early succession is much more unpredictable as most strategies can establish on the bare substrate (see Fig. 9 in Appendix). Plant biomass increases along succession while both soil P and N availabilities decrease. This competition-driven scenario is close to what [START_REF] Connell | Mechanisms of succession in natural communities and their role in community stability and organization[END_REF] called the tolerance model.

Nutrient supplies not only control succession qualitatively (competition-vs facilitation-driven), they also influence it quantitatively. This can be read on Fig. 5 using the envelope and the impact vectors from the resource-ratio theory. First, because every succession endpoint resource level is located on the envelope, final soil N and P are negatively correlated.

For example, the higher the N supply, the higher the N and the lower the P levels at the succession endpoint. The system responds differently to P substrate availability depending on the succession scenario. When competition-driven, an increase in S P leads to an increase in final soil P; the opposite happens in the facilitation case, as a P-richer substrate eventually leads to a more P-depleted soil. These counter-intuitive results stem from the community-level regulation of N availability around N * c , a consequence of the coupling between the ecological positive feedback loop and community assembly. When taken separately, each plant species in the facilitation-driven scenario is only P-limited, as an increase in S N does not lead to an increase in biomass. This contrasts with the community-wide response to nutrient addition, as the community endpoint biomass responds positively to an increase in both S N and S P , due to species turnover.

Discussion

In this paper, we show how facilitation-driven succession can emerge from repeated invasions of a diversity of N-fixing species in interaction with soil N and P. We characterize the traits, nutrient and biomass dynamics of the developing ecosystem during facilitation-driven succession and show that late succession presents an increasing sensitivity to catastrophic shifts. Finally, we show how changes in N and P external supplies influence succession characteristics and even succession scenario, shifting from facilitation-to competition-driven.

Facilitation-driven succession differs from classic competition-driven succession by the fact that the invasion success of non-pioneer species relies on the successful invasion and environmental modifications of previous strategies. Such a mechanism gives a straightforward way to test for facilitation-driven succession empirically, as initially proposed by [START_REF] Connell | Mechanisms of succession in natural communities and their role in community stability and organization[END_REF]: experimentally remove the pioneer species and watch for changes in the succession dynamics.

This approach gave evidence for facilitation-driven succession in rocky intertidal communities [START_REF] Farrell | Models and mechanisms of succession: an example from a rocky intertidal community[END_REF], but seemed to rule out this mechanism on old-field secondary succession [START_REF] Hils | Species removals from a first-year old-field plant community[END_REF][START_REF] Armesto | Removal experiments to test mechanisms of plant succession in oldfields[END_REF] and in recent experiments on carrion [START_REF] Michaud | Facilitation may not be an adequate mechanism of community succession on carrion[END_REF].

Our model however suggests that such tests should not be restricted to pioneer species, as any species can facilitate the invasion of subsequent strategies. Every species can thus play a role during facilitation-driven succession and slow it down when absent. Indeed, the mean succession time of facilitation-driven succession is expected to depend on the colonization rates of all the species in the pool, in stark contrast with competition-driven succession where the mean succession time is essentially the inverse of the colonization rate of the climax strategy, usually much shorter. We therefore suggest that competition versus facilitation-driven succession could be tested experimentally by rarefying the species pool: a subsequent increase in mean succession time would support facilitation-driven succession. If rarefying opens a large gap in the species pool, succession could get blocked in its early stages. However, redundancy between closely functioning N-fixers ensures a certain robustness of facilitation-driven succession to species loss. This illustrates a dynamical 'insurance effect' on ecosystem development provided by the diversity of N-fixing strategies of the regional species pool [START_REF] Yachi | Biodiversity and ecosystem productivity in a fluctuating environment: The insurance hypothesis[END_REF].

Our model describes a replacement from high-fixers to low-fixers during succession, consistently with some empirical evidence reviewed by [START_REF] Vitousek | Biological nitrogen fixation: rates, patterns and ecological controls in terrestrial ecosystems[END_REF]. This is accompanied with an increasing sensitivity to catastrophic shifts of facilitation-driven succession as it approaches its climax, inevitably leading to a late succession collapse towards the bare substrate in the presence of disturbances. From there, succession has to start again from scratch, meaning that in practice, facilitation-driven succession exhibits periodic alternations between developmental phases and sudden collapses, like [START_REF] Watt | Pattern and process in the plant community[END_REF] autogenic succession cycles or other community assembly cycles [START_REF] Morton | Regional species pools and the assembly of local ecological communities[END_REF][START_REF] Steiner | Cyclic assembly trajectories and the generation of scale-dependent productivitybiodiversity relationships[END_REF]). Such cycles have been observed on empirical systems such as individual associations between cacti and nursing shrubs [START_REF] Yeaton | The dynamics of a succulent karoo vegetation[END_REF][START_REF] Valiente-Banuet | Interaction between the cactus Neobuxbaumia tetetzo and the nurse shrub Mimosa luisana[END_REF]), but we are not aware of any empirical evidence at the ecosystem scale.

From a theoretical perspective, such increasing sensibility to collapse is a consequence of the opposition between the private costs of facilitation and their collective benefits, which selects for cheaters with less efficient fixation, an ecological phenomenon known as the tragedy of the commons [START_REF] Hardin | The tragedy of the commons[END_REF]) or evolutionary suicide [START_REF] Ferrière | Adaptive responses to environmental threats: evolutionary suicide, insurance, and rescue[END_REF][START_REF] Gyllenberg | Necessary and sufficient conditions for evolutionary suicide[END_REF]. Some mechanisms such as spatial aggregation [START_REF] Kéfi | Evolution of local facilitation in arid ecosystems[END_REF], kin selection [START_REF] Lion | Self-structuring in spatial evolutionary ecology[END_REF] or 'leakiness' from the 'Black Queen Hypothesis' [START_REF] Morris | Black Queen evolution: the role of leakiness in structuring microbial communities[END_REF] could explain how facilitation is maintained during succession, potentially reconciling our approach with [START_REF] Odum | The strategy of ecosystem development[END_REF] hypothesis of maximal stability at the climax. In fact, some ecosystems such as tropical forests present an abundance of N-fixing plants despite high availability of N [START_REF] Hedin | The nitrogen paradox in tropical forest ecosystems[END_REF]. Accounting for mechanisms such as the ones mentioned above is thus necessary to understand ecosystem stability and the maintenance of facilitation in late succession.

Another prediction of our model is that the succession scenario depends on the external nutrient supplies. This can be tested through nitrogen-addition experiments, by characterizing succession under various levels of fertilization of an N-depleted substrate. We expect facilitation-driven succession to vanish under high levels of fertilization, and mean succession time and final soil P concentration to decrease along fertility gradients. Conversely, the absence of such patterns would suggest that the initial substrate is too rich for facilitation to drive the succession, as exemplified by some studies of secondary succession in old-fields [START_REF] Huberty | Effects of nitrogen addition on successional dynamics and species diversity in Michigan old-fields[END_REF][START_REF] Wilson | Quadratic variation in old-field species richness along gradients of disturbance and nitrogen[END_REF]. In this perspective, we predict that the recent global increase in N availability due to human activities (Vitousek et al. 1997) has favored competition-driven succession and decreased mean succession times.

In our approach, we parsimoniously assumed equal colonization probability for every strategy. However, species differ in their dispersal and competitive abilities, e.g. due to seed size/number trade-offs [START_REF] Turnbull | Seed mass and the competition/colonization trade-off: a sowing experiment[END_REF], which leads to competition/colonization trade-offs at the core of the corresponding succession scenario [START_REF] Tilman | Competition and biodiversity in spatially structured habitats[END_REF]. How would a classical competition/colonization trade-off combine with the facilitation-driven succession scenario? This would depend on how maximal fixation rate correlates with dispersal ability. General information on these correlations is lacking, but N-fixing pioneers tend to be very good dispersers, e.g. microorganisms that travel worldwide attached to dust particules [START_REF] Gorbushina | Life in Darwin's dust: intercontinental transport and survival of microbes in the nineteenth century[END_REF] or Dryas drummondii and its wind-dispersed plumed achenes [START_REF] Lawrence | The role of Dryas drummondii in vegetation development following ice recession at Glacier Bay, Alaska, with special reference to its nitrogen fixation by root nodules[END_REF]. Such dispersal syndromes could in turn explain how these inferior competitors -crucial to initiate succession -are maintained at the landscape level in a metacommunity perspective [START_REF] Hastings | Disturbance, coexistence, history, and competition for space[END_REF][START_REF] Leibold | The metacommunity concept: a framework for multi-scale community ecology[END_REF]).

Our approach relied on two critical assumptions, justified by previous works, that fixation is facultative and that it incurs a cost in terms of P acquisition. If we instead considered obligate N-fixers, we expect facilitation through enriched soil nitrogen to be much stronger. This would dramatically accelerate facilitation-driven succession, and allow most Nfixing species to invade after a short pioneer-dominated phase. This, in turn, could further weaken late succession, and lead to an almost certain collapse of the ecosystem when low-fixing strategies, too weak to maintain the positive feedback loop, happen to invade. This, again, emphasizes the need for additional mechanisms to understand late succession stability.

Loosening the trade-off between that maximal N-fixation and P-acquisition does not qualitatively alter our analysis (Appendix D). Completely relaxing the trade-off would, in the current version of the model, turn the total-fixer into a 'Hutchinsonian demon' [START_REF] Kneitel | Trade-offs in community ecology: linking spatial scales and species coexistence[END_REF]), a species that would exclude all the others and freeze succession.

The trade-off we use relies on the idea of fixed allocation between two plant traits only. Yet, it is possible that other limiting factors involved in succession actually limit fixation rates. For example, herbivory is known to affect N-fixers and succession [START_REF] Ritchie | Responses of legumes to herbivores and nutrients during succession on a nitrogenpoor soil[END_REF][START_REF] Kurokawa | Plant traits, leaf palatability and litter decomposability for co-occurring woody species differing in invasion status and nitrogen fixation ability[END_REF]. Adding a third trait to our approach would lead to a multidimensional problem, providing a promising direction for future modeling work on facilitation-driven succession.

We focused on N-mediated facilitation, but there are other ways through which early colonists can facilitate later invading species [START_REF] Bruno | Inclusion of facilitation into ecological theory[END_REF]), e.g. by improving water infiltration [START_REF] Klausmeier | Regular and irregular patterns in semiarid vegetation[END_REF][START_REF] Hillerislambers | Vegetation pattern formation in semi-arid grazing systems[END_REF], providing shade to light-intolerant species [START_REF] Gerla | Photoinhibition and the assembly of light-limited phytoplankton communities[END_REF], or increasing nutrient availability through improved soil retention or changes in soil pH [START_REF] Hinsinger | P for two, sharing a scarce resource: soil phosphorus acquisition in the rhizosphere of intercropped species[END_REF]. Such mechanisms could explain decreased nutrient leaching rates during succession [START_REF] Odum | The strategy of ecosystem development[END_REF][START_REF] Vitousek | Ecosystem succession and nutrient retention: a hypothesis[END_REF] and the exclusion of N-fixing strategies in some late successions [START_REF] Chapin | Mechanisms of primary succession following deglaciation at Glacier Bay, Alaska[END_REF][START_REF] Menge | Evolutionary tradeoffs can select against nitrogen fixation and thereby maintain nitrogen limitation[END_REF], two situations that cannot happen within our model.

To conclude, we presented two principal results regarding facilitation-driven succession by N-fixers. First, we showed that resource availability determines if the succession scenario is facilitation-or competition-driven. Second, we showed that late facilitation-driven succession is sensitive to catastrophic shifts, leading to succession cycles. We therefore expect that the addition of some mechanisms, yet to be explored, is crucial to explain both the maintenance of facilitation during late succession and the stability of these ecosystems. (dark blue) and maximal N-fixation (light blue) combine in Liebig's law to determine the growth g tot , uptake and fixation rates (black arrow) of facultative N-fixers along a N availability gradient. When soil P is more limiting than soil N, plants are P-limited and fixation is unnecessary (far left). When soil N is more limiting than soil P, there are three possibilities: either maximal fixation rate is large enough for N not to be limiting, leading to P-limitation with partial N-fixation (center left), or maximal fixation is just right to reach colimitation (center right), or fixation is insufficient to overcome N limitation, leading to N-limitation with maximal N-fixation (far right). The mathematical expression of the effective fixation rate associated with these different scenarios can be found in Appendix A. The ZNGI corners all line up on a decreasing envelope (dashed grey), which materializes the trade-off between P * and N * nf . (b),(c) and (d): Ecological phase diagrams along N and P availabilities for three different N-fixing strategies: low, high and total N-fixer corresponding respectively to F/m = 0.15, 0.6 and 1. The ZNGI (thick, L-shaped), a strategy's minimal nutrient requirements for growth, also gives the possible nutrient levels at equilibrium, which can be either stable (blue, plain) or unstable (purple, dashed). The impact vectors (blue and purple vectors) correspond to a strategy's per capita impact on the two nutrients. They point down when the impact on N is negative (consumption), up when it is positive (accumulation). The blue-shaded region (labeled B F ) gives the region of external nutrient supplies S = (S P , S N ) (e.g. black dot) under which the N-fixer can subsist, leading to a vegetated ecosystem. In the white region (labeled ∅), supplies are too low and the N-fixer goes extinct, leading to a bare substrate. The purple-shaded region (labeled B F |∅) has the vegetated ecosystem and the bare substrate as two alternative stable states. This can be seen in the inset in (c), where the biomass at equilibrium (blue) of these two stable states (plain) and the unstable state that separates them (dashed) have been plotted against a nutrient gradient with increasing S N and fixed S P (black, vertical line). Biomass is constant along this N gradient because P is limiting ans S P kept constant. The ZNGIs of two high N-fixing strategies F 1 and F 2 are superimposed, following the color code of Fig2. F 1 is a better N-fixer than F 2 (F 1 > F 2 ), but worse P-competitor. Only F 1 can invade the bare substrate (S), the supply point (black dot) being located above the invasion N requirements of F 1 (dashed, purple) but below the ones of F 2 . As F 1 establishes, it moves the nutrient levels to R 1 , where invasion by F 2 is possible. As F 2 establishes, it further moves the nutrient levels to R 2 , competitively excluding F 1 . As S belongs in F 2 alternative stable state region (purple), F 2 can subsist there even after displacing F 1 . (b) Same successional step represented on a PIP, picturing the invasion success of an invader (F ) versus a resident state (F and ∅). The resident state can either be the bare substrate (∅), or a previously established strategy with trait F that is globally stable (white, e.g. F 1 ) or sensitive to catastrophic shift because the bare substrate is an alternative stable state (gray, e.g. F 2 ). Some strategies F are unable to be resident (black). Positive invasion pairs, by an invader F in a resident F or ∅, are pictured in blue. F 1 can invade the bare substrate while F 2 cannot. After F 1 establishes, F 2 can invade F 1 and displace it, leading to a bistable resident state. All the trajectories were the result of 100 iterations of random colonization attempts from the species pool with equal probability followed by an ecological sorting phase resulting in either 1) extinction of the invader; 2) replacement of the resident by the invader. (a) and (b): same situation as Fig. 3 iterated 100 times. Rapidly, as the trajectory gets close to the ESS, virtually every colonization attempt fails, which explains why only 5 colonizations succeed out of 100. (c), (d), (e) and (f): ecosystem properties dynamics along succession, respectively maximal fixation rate F , biomass BF and soil nutrient levels P and N . 10 random trajectories are represented (colored, plain) as well as the three quartiles (thick black, median plain; 1 st and 3 rd quartiles dashed) of the full distribution of probability of the succession process, obtained with a Markov chain. F and P decrease during succession, while BF and N increase, all of them reaching a climax as we get close to the ESS. Contrary to competition-driven succession (Fig. 9 in Appendix E), facilitation ensures that ecosystem properties are 'channeled' during early succession.

= l N (S N -N ) -q P [g N ( N ) -λm] BF (11)
There is no explicit analytical solution for a general growth function g N . However, the solution can still be drawn graphically using the supply point map as was done in the main text.

Substituting these results in the condition for P-limitation with N-fixation g N ( N ) < g P ( P ) < g N ( N ) + F , we rewrite

this condition N * nf > N > N * f with N * f = g -1 N (m-F ).
To get a more explicit expression of this condition as a function of the parameters of the system, we introduce the parameter β = (1-λ)m/[(1-λ)m-F ]. When F < (1-λ)m, we call the plant a low N-fixer and β > 0. The condition for P-limitation with N-fixation then rewrites βγ(S N -N * f ) > S P -P * > γ(S N -N * nf ).

When F > (1 -λ)m, we call the plant a high N-fixer and β < 0. The condition for P-limitation with N-fixation then

rewrites S P -P * > max[γ(S N -N * nf ), βγ(S N -N * f )].
A last condition for this equilibrium to exist is BF > 0, which always translates into S P -P * > 0.

Also note that we defined in the main text the threshold N * c in nitrogen concentration under which the N-fixer has a positive impact on soil nitrogen. Formally, it satisfies g N (N * c ) -λm = 0, which gives N * c = g -1 N (λm).

N-limited N-fixing plant: This case corresponds to the far right panel of Fig. 1b. Formally, the N-limitation with N-fixation for B corresponds to the constraint g N ( N ) + F < g P ( P ). We get from eq. ( 2a) that N = N * f . BF is then deduced from eq. ( 2b) and is equal to:

BF = l N q N S N -N * f (1 -λ)m -F (12)
and from eq. (2c) we get:

P = S P - q P l P (1 -λ)m BF (13)
Substituting these results in the condition for N-limitation with N-fixation g N ( N ) + F < g P ( P ), we rewrite this condition P > P * or more explicitly S P -P * > βγ(S N -N * f ). A last condition for this equilibrium to exist is BF > 0, which translates into βγ(S N -N * f ) > 0, i.e. S N -N * f > 0 for low N-fixers and S N -N * f < 0 for high N-fixers.

These results can be summarized in the supply plane by counting the different equilibria associated with each supply point (S P , S N ). In the case of a low N-fixer (β > 0), the P-limited non-fixing, P-limited N-fixing and N-limited N-fixing supply regions do not overlap and form together the region of plant presence -defined by the overall conditions S P > P * and S N > N * f . In the case of a high N-fixer (β < 0), the N-limited N-fixing region is included in the P-limited N-fixing region, which means that these supply points are associated with two possible equilibria with the plant present. We show in the next section that this leads to alternative stable states.

Stability

The Jacobian of eq. ( 2) writes:

J (B, P, N ) =       g tot -m ∂ P g tot B ∂ N g tot B i P -l P + ∂ P i P B ∂ N i P B i N ∂ P i N B -l N + ∂ N i N B       (14)
where the notation i P (P, N ) = -q P [g tot (P, N ) -λm] and i N (P, N ) = -q N [min[g P (P ), g N (N )] -λm] designate the per capita net impacts of N-fixers on P and N.

Bare substrate ∅: With BF = 0, N = S N and P = S P , we have:

J ( BF , P , N ) =       min[g P (S P ), g N (S N ) + F ] -m 0 0 i P -l P 0 i N 0 -l N       (15)
the three eigenvalues are located on the diagonal; we conclude from the condition min[g P (S P ), g N (S N ) + F ] -m < 0 that the bare substrate ∅ is locally stable only when the plant cannot invade, i.e. for supply points located under the L-shaped ZNGI.

When the plant is present, the Jacobian J rewrites:

J ( BF , P , N ) =       0 ∂ P g tot BF ∂ N g tot BF i P -l P + ∂ P i P BF ∂ N i P BF i N ∂ P i N BF -l N + ∂ N i N BF       (16)
P-limited non-fixing plant: this leads to:

J ( BF , P , N ) =       0 ∂ P g P BF 0 -q P (1 -λ)m -l P -q P ∂ P g P BF 0 -q N (1 -λ)m -q N ∂ P g P BF -l N       (17)
which means than λ N = -l N , λ -λ + = q P (1 -λ)m∂ P g P BF > 0 and λ -+ λ + = -l P -q P ∂ P g P BF < 0; thus, all the eigenvalues are negative and this equilibrium is always locally stable.

P-limited N-fixing plant: this leads to:

J ( BF , P , N ) =       0 ∂ P g P BF 0 -q P (1 -λ)m -l P -q P ∂ P g P BF 0 -q N (g N -λm) 0 -l N -q N ∂ N g N BF       (18)
which, pretty similarly to the previous case, means than λ N = -l N -q N ∂ N g N B, λ -λ + = q P (1 -λ)m∂ P g P B > 0 and λ -+ λ + = -l P -q P ∂ P g P B < 0; thus, all the eigenvalues are negative and this equilibrium is always locally stable.

N-limited N-fixing plant: this leads to:

J ( BF , P , N ) =       0 0 ∂ N g N BF -q P (1 -λ)m -l P -q P ∂ N g N BF -q N (g N -λm) 0 -l N -q N ∂ N g N BF       (19)
which means that λ

P = -l P , λ -λ + = q N [(1 -λ)m -F ]∂ N g N BF and λ -+ λ + = -l N -q N ∂ N g N BF < 0.
The sign of λ -λ + is positive if and only if β < 0, i.e. for high N-fixers, which means that the N-limited equilibria are always unstable in this case.

To conclude, the P-limited equilibria are always stable when they exist. The N-limited equilibria are stable for low N-fixers (β > 0), and unstable for high N-fixers (β < 0). The bare substrate equilibria are stable under the ZNGI, i.e.

whenever S P < P * or S N < N * f . In the case of a low N-fixer, the regions of stable bare substrate and stable plant presence are complementary so that there is a unique stable equilibrium associated with each supply point. In the case of a high N-fixer, the N-limited N-fixing plant region located under the ZNGI -formally defined by S N < N * f and S P -P * > βγ(S N -N * f ) -corresponds to alternative stable states, as the unstable N-limited N-fixing equilibrium separates a stable bare substrate equilibrium from a stable P-limited N-fixing equilibrium.

parameters, and the associated units. We also give the specific numerical values and functional forms used to generate the figures, and recapitulate the particular expression of the main lumped parameters. rate F at the ESS compared to the linear case. Looking at the trajectories themselves, one can see that maximal fixation rates indeed decrease faster compared to the linear trade-off case (Fig. 8c). The changes in ecosystem properties during succession happen even faster (Fig. 8d-f). This comes from the fact that the P-acquisition of mid-succession strategies with intermediate F are already very close to the ones of the ESS as, a consequence of the weak trade-off. To conclude, the 650 main effect of a looser trade-off between N-fixation and P-acquisition is to accelerate the rate of succession. We however note that, next to these quantitative changes in the results, the global dynamics remain qualitatively analogous to the one detailed in the main text. e) and (f): ecosystem properties dynamics along succession, respectively maximal fixation rate F , biomass BF and soil nutrient levels P and N . 10 random trajectories are represented (colored, plain) as well as the three quartiles (thick black, median plain; 1 st and 3 rd quartiles dashed) of the full distribution of probability of the succession process, obtained with a Markov chain. F , P and N decrease during succession while BF increases, all of them reaching a climax as we get close to the ESS. Contrary to facilitation-driven succession (Fig. 4), the majority of the N-fixing species can establish in the bare substrate.

Figure 1 :

 1 Figure 1: (a) Flow diagram of the model. (b) Schematic representation of how soil P availability (orange), soil N availability (dark blue) and maximal N-fixation (light blue) combine in Liebig's law to determine the growth g tot , uptake and fixation rates (black arrow) of facultative N-fixers along a N availability gradient. When soil P is more limiting than soil N, plants are P-limited and fixation is unnecessary (far left). When soil N is more limiting than soil P, there are three possibilities: either maximal fixation rate is large enough for N not to be limiting, leading to P-limitation with partial N-fixation (center left), or maximal fixation is just right to reach colimitation (center right), or fixation is insufficient to overcome N limitation, leading to N-limitation with maximal N-fixation (far right). The mathematical expression of the effective fixation rate associated with these different scenarios can be found in Appendix A.

F‰Figure 2 :

 2 Figure 2: (a) set of ZNGIs (blue, L-shaped) representing the competitive abilities of N-fixers from the species pool with fixation abilities F/m comprised between 0 and 1. The ZNGI corners all line up on a decreasing envelope (dashed grey), which materializes the trade-off between P * and N * nf . (b),(c) and (d): Ecological phase diagrams along N and P availabilities for three different N-fixing strategies: low, high and total N-fixer corresponding respectively to F/m = 0.15, 0.6 and 1. The ZNGI (thick, L-shaped), a strategy's minimal nutrient requirements for growth, also gives the possible nutrient levels at equilibrium, which can be either stable (blue, plain) or unstable (purple, dashed). The impact vectors (blue and purple vectors) correspond to a strategy's per capita impact on the two nutrients. They point down when the impact on N is negative (consumption), up when it is positive (accumulation). The blue-shaded region (labeled B F ) gives the region of external nutrient supplies S = (S P , S N ) (e.g. black dot) under which the N-fixer can subsist, leading to a vegetated ecosystem. In the white region (labeled ∅), supplies are too low and the N-fixer goes extinct, leading to a bare substrate. The purple-shaded region (labeled B F |∅) has the vegetated ecosystem and the bare substrate as two alternative stable states. This can be seen in the inset in (c), where the biomass at equilibrium (blue) of these two stable states (plain) and the unstable state that separates them (dashed) have been plotted against a nutrient gradient with increasing S N and fixed S P (black, vertical line). Biomass is constant along this N gradient because P is limiting ans S P kept constant.

Figure 3 :

 3 Figure 3: Representation of a facilitation-driven successional step in the nutrient (a) and strategy space (b), emerging under low N-and high P-availability, i.e. S = (S P , S N ) = (1.8, 0.05). (a)The ZNGIs of two high N-fixing strategies F 1 and F 2 are superimposed, following the color code of Fig2. F 1 is a better N-fixer than F 2 (F 1 > F 2 ), but worse P-competitor. Only F 1 can invade the bare substrate (S), the supply point (black dot) being located above the invasion N requirements of F 1 (dashed, purple) but below the ones of F 2 . As F 1 establishes, it moves the nutrient levels to R 1 , where invasion by F 2 is possible. As F 2 establishes, it further moves the nutrient levels to R 2 , competitively excluding F 1 . As S belongs in F 2 alternative stable state region (purple), F 2 can subsist there even after displacing F 1 . (b) Same successional step represented on a PIP, picturing the invasion success of an invader (F ) versus a resident state (F and ∅). The resident state can either be the bare substrate (∅), or a previously established strategy with trait F that is globally stable (white, e.g. F 1 ) or sensitive to catastrophic shift because the bare substrate is an alternative stable state (gray, e.g. F 2 ). Some strategies F are unable to be resident (black). Positive invasion pairs, by an invader F in a resident F or ∅, are pictured in blue. F 1 can invade the bare substrate while F 2 cannot. After F 1 establishes, F 2 can invade F 1 and displace it, leading to a bistable resident state.

Figure 4 :

 4 Figure4: Facilitation-driven succession trajectories, under low N and high P supplies, i.e. S = (S P , S N ) = (1.8, 0.05). All the trajectories were the result of 100 iterations of random colonization attempts from the species pool with equal probability followed by an ecological sorting phase resulting in either 1) extinction of the invader; 2) replacement of the resident by the invader. (a) and (b): same situation as Fig.3iterated 100 times. Rapidly, as the trajectory gets close to the ESS, virtually every colonization attempt fails, which explains why only 5 colonizations succeed out of 100. (c), (d), (e) and (f): ecosystem properties dynamics along succession, respectively maximal fixation rate F , biomass BF and soil nutrient levels P and N . 10 random trajectories are represented (colored, plain) as well as the three quartiles (thick black, median plain; 1 st and 3 rd quartiles dashed) of the full distribution of probability of the succession process, obtained with a Markov chain. F and P decrease during succession, while BF and N increase, all of them reaching a climax as we get close to the ESS. Contrary to competition-driven succession (Fig.9in Appendix E), facilitation ensures that ecosystem properties are 'channeled' during early succession.

Figure 5 :FFigure 6 :

 56 Figure 5: Community-level phase diagram summarizing the different succession scenarios along external nutrient supplies S P and S N . When N and P supplies are too low (white region on the left of the gray curve), not a single strategy can subsist and the ecosystem stays in the bare substrate state (0). Low N supply (S N < N * c ) favors high N-fixation rates with positive net impact on N, which leads to facilitation-driven succession (I). High N supply (S N > N * c ) favors low N-fixation rate that do not compensate for N uptake, which leads to competition-driven succession (III).

  0

Figure 8 :Figure 9 :

 89 Figure 8: Facilitation-driven succession trajectories in the case of a weak trade-off between N-fixation and P-acquisition. The rest of the situation is strictly similar to Fig. 4.

Table 1 :

 1 Model variables, functions and parameters Trade-off function α P (F ) = α P,0 + h(F )(α P,1 -α P,0 )

	Symbol	Meaning	Value used	Units
	Dynamical variables
	P	Phosphorous (P) availability	-	gP . m -2
	N	Nitrogen (N) availability	-	gN . m -2
	B F	N-fixer density	-	g of biomass . m -2
	Ecological parameters
	F	N-fixer maximal fixation rate	0 -0.5	yr -1
	m	N-fixer per capita death rate	1	yr -1
	α P	N-fixer P acquisition rate	see below	gP -1 . m 2 . yr -1
	α N	N-fixer N acquisition rate	2	gN -1 . m 2 . yr -1
	q P	N-fixer P-to-biomass ratio	1	gP . (g of biomass) -1
	q N	N-fixer N-to-biomass ratio	1	gN . (g of biomass) -1
	I P	P external input rate	0 -2	gP . m -2 . yr -1
	I N	N external input rate	0 -0.8	gN . m -2 . yr -1
	l P	P leaching rate	1	yr -1
	l N	N leaching rate	3	yr -1
	λ	Recycling efficiency	0.75	-
	Growth functions
	g P (P ) = α P P	N-fixer growth rate when soil P is limiting	-	yr -1
	g N (N ) = α N N	N-fixer growth rate when soil N is limiting	-	yr -1
		P acquisition rate	-	gP -1 . m 2 . yr -1
	α P,0	Minimal P acquisition rate	0.666	gP -1 . m 2 . yr -1
	α P,1	Maximal P acquisition rate	2.666	gP -1 . m 2 . yr -1
	h(F ) = 1 -F/m	linear trade-off function	-	-
	Lumped parameters
	S P = I P /l P	P supply	0 -2	gP . m -2
	S N = I N /l N	N supply	0 -0.8	gN . m -2
	P * = m/α P	N-fixer minimal P requirements	-	gP . m -2
	N * nf = m/α N N * f = (m -F )/α N N * c = λm/α N	N-fixer minimal N requirements if not fixing N-fixer minimal N requirements when fixing N-fixer positive impact threshold	0.5 -0.375	gN . m -2 gN . m -2 gN . m -2

Acknowledgments

Vanessa Ezenwa and two anonymous reviewers provided helpful comments on the manuscript. We thank Bertrand Cloez, Christopher Klausmeier, Yannick Kouakep Tchaptchié, Elena Litchman and the other members of the Litchman-Klausmeier Lab for useful discussions on this subject. This is Kellogg Biological Station contribution #2011.

A Effective fixation rate

In this appendix, we show how the modeling assumptions introduced in the main text on facultative fixation lead to a downregulation of the effective fixation rate, denoted φ, by N and P availabilities. Doing so, we present an alternativeand equivalent -formulation of the plant-soil model that is more explicit towards the amount of nitrogen that is actually fixed by the plant.

As a starting point, we recall the two dynamical equations describing N-fixing population biomass B F and soil nitrogen concentration N :

The idea here is to remember that we assumed both mass balance and fixed stoichiometry q N . Because of this, the effective fixation rate φ -the rate at which N is actually fixed -is simply equal to the difference between the plant growth rate and the plant N uptake rate, the latter divided by the stoichiometric coefficient q N :

From there, we can use the mathematical identity min(a, b + c) -min(a, b) = max[0, min(a -b, c)] to obtain the more compact expression:

This effective fixation rate φ is represented in Fig. 6 as a function of soil N limitation g P (P ) -g N (N ), and quantifies the fixation mechanism described in Fig. 1b. When P is limiting (g P (P ) -g N (N ) < 0), there is no need to fix and φ = 0. When N is limiting while the plant fixes at its maximal rate (g P (P ) -g N (N ) > F ), we logically have φ = F .

In the region in-between (0 < g P (P ) -g N (N ) < F ), the plant fixes just the amount needed to reach co-limitation, i.e. φ = g P (P ) -g N (N ).

To conclude, we have obtained an alternative -and equivalent -way to write the model of eq. ( 2):

B Analytical results of the plant-soil model

In this appendix, we present the analytical calculations underlying the graphical results of the plant-soil model displayed in the main text. More precisely, we compute the different equilibria of the plant-soil model and analyze the conditions for their stability as functions of phosphorous and nitrogen supplies.

The results presented here hold for general growth functions g X (X) where X = P, N given that they satisfy ∂g X /∂X and lim X→∞ g X (X) > m.

Equilibria

There are four cases to consider : the bare substrate, which consists in the trivial equilibrium in the absence of plant, and three situations associated with the presence of plants -P-limited non-fixing, P-limited N-fixing and N-limited N-fixing (see Fig. 1b).

Bare substrate ∅: This case corresponds to BF = 0. From eq. ( 2a), it follows that N = S N = I N /l N and P = S P = I P /l P . This equilibrium always exists; we will see in the next subsection the conditions for its stability.

P-limited non-fixing plant: This case corresponds to the far left panel of Fig. 1b. Formally, the P-limitation without N-fixation for B F corresponds to the constraint g P ( P ) < g N ( N ) on the nutrient levels at equilibrium. We will identify later the parameters leading to this condition. Eq. (2a) at equilibrium then gives P = P * = g -1 P (m). BF is then deduced from eq. ( 2b) and is equal to:

and from eq. (2c) we get:

Substituting these results in the condition for P-limitation without N-fixation g P ( P ) < g N ( N ), this condition rewrites

, or more explicitly γ(S N -N * nf ) > S P -P * with γ = l N q P /l P /q N . A last condition for this equilibrium to exist is BF > 0, which simply translates into S P -P * > 0.

P-limited N-fixing plant: This case corresponds to the center left panel of Fig. 1b. Formally, the P-limitation with Nfixation for B corresponds to the constraint g N ( N ) < g P ( P ) < g N ( N ) + F on the nutrient levels at equilibrium. Similarly to the previous case, we get from eq. ( 2a) that P = P * = g -1 P (m) and BF is equal to:

C Model variables, functions and parameters

In this appendix, we recapitulate the notation used throughout the paper in terms of dynamical variables and ecological D Weak trade-off between N-fixation and P-acquisition

In this appendix, we explore the effect of loosening the trade-off between N-fixation and P-acquisition introduced in the main text. This is achieved by increasing the concavity of the relationship between N-fixation and P-acquisition. We show that such a looser trade-off results in an acceleration of facilitation-driven succession and a slight decrease in the maximal fixation rate of the ESS.

The figures of the main text were generated using a linear trade-off between N-fixation and P-acquisition, of the form α P (F ) = α P,0 + h(F )(α P,1 -α P,0 ), with linear h(F ) = 1 -F/m (see Table 1). We now assume that h is instead a concave function of the form:

where the parameter quantifies the concavity of the trade-off. When = 1, we recover the linear trade-off (blue line, Fig. 7), while < 1 (taken equal to 1/10 in this appendix) gives a concave trade-off (purple line, Fig. 7). We can see why such a concave relationship leads to a weaker trade-off between N-fixation and P-acquisition: 1) because the concave trade-off is located above the linear trade-off, the strategies located on the concave trade-off are overall more competitive than the ones located on the linear trade-off; 2) the strategies located at the top-right "corner" of the trade-off are able to achieve both high N-fixation and high P-acquisition rates.

We generated facilitation-driven succession trajectories similar to Fig. 4 in the case of this concave trade-off. The results are summarized in Fig. 8. Several points can be highlighted. As expected, the concave trade-off leads to an increased convexity of the ZNGI envelope, graphically translating the fact that some strategies are good at both N-fixation and P-acquisition (Fig. 8a). On the PIP, we see that the positive invasion region is broadened, which accelerates further