Gaël Bardon
email: gael.bardon2@gmail.com

Frédéric Barraquand
email: frederic.barraquand@u-bordeaux.fr

Effects of stage structure on coexistence: mixed benefits

Keywords: coexistence, stage structure, competition, matrix population models

The properties of competition models where all individuals are identical are relatively well-understood; however, juveniles and adults can experience or generate competition differently. We study here less wellknown structured competition models in discrete time that allow multiple life history parameters to depend on adult or juvenile population densities. A numerical study with Ricker density-dependence suggested that when competition coefficients acting on juvenile survival and fertility reflect opposite competitive hierarchies, stage structure could foster coexistence. We revisit and expand those results.

First, through a Beverton-Holt two-species juvenile-adult model, we confirm that these findings do not depend on the specifics of density-dependence or life cycles, and obtain analytical expressions explaining how this coexistence emerging from stage structure can occur. Second, we show using a community-level sensitivity analysis that such emergent coexistence is robust to perturbations of parameter values. Finally, we ask whether these results extend from two to many species, using simulations. We show that they do not, as coexistence emerging from stage structure is only seen for very similar life-history parameters. Such emergent coexistence is therefore not likely to be a key mechanism of coexistence in very diverse ecosystems, although it may contribute to explaining coexistence of certain pairs of intensely competing species.

Introduction

A basic tenet of demography and population ecology is that species vital rates, such as survival and fertility, can vary with age or stage [START_REF] Leslie | On the use of matrices in certain population mathematics[END_REF][START_REF] Cushing | An introduction to structured population dynamics[END_REF][START_REF] Caswell | Matrix Population Models: Construction, Analysis, and Interpretation[END_REF]. Surprisingly, little of this rich agestructured population-level theory influences community-level coexistence theory, that by contrast mostly builds on unstructured Lotka-Volterra competition or consumer-resource models [START_REF] Chesson | Mechanisms of maintenance of species diversity[END_REF][START_REF] Barabás | The effect of intra-and interspecific competition on coexistence in multispecies communities[END_REF][START_REF] Letten | Linking modern coexistence theory and contemporary niche theory[END_REF]. There have been, however, some calls to include more demography into community ecology (e.g., [START_REF] Miller | Thinking inside the box: community-level consequences of stagestructured populations[END_REF] and coexistence studies in particular, as well as a number of empirically-driven modelling studies doing so [START_REF] Péron | Integrated modeling of communities: parasitism, competition, and demographic synchrony in sympatric ducks[END_REF][START_REF] Chu | Large niche differences emerge at the recruitment stage to stabilize grassland coexistence[END_REF]. They show in general that the vital rates of the various life-stages do not react in the same way to changes in the densities of juvenile and adult life-stages. These differences may well create new avenues for coexistence.

Although mainstream coexistence theory often neglects the complexity of life histories, some stagestructured coexistence theory has been previously developed: several theoretical models, in past decades [START_REF] Haigh | Can there be more predators than prey?[END_REF][START_REF] Loreau | Competitive exclusion and coexistence of species with complex life cycles[END_REF] as well as in more recent literature [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF][START_REF] Fujiwara | Coexistence of competing stage-structured populations[END_REF] have considered complex life cycles in which ontogenetic changes occur, leading to models where vital rates can be differentially affected by competition. Such complexity can promote equalizing mechanisms [START_REF] Fujiwara | Coexistence of competing stage-structured populations[END_REF] as well as niche differentiation [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF]. A main observation of [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF] was the occurrence of so-called 'emergent coexistence' in two-species two-stages models, which was defined as a coexistence equilibrium where competition on a single vital rate would lead to exclusion, but the combination of competition processes affecting two different vital rates leads to coexistence (e.g., competition affecting juvenile survival leads to species 1 winning, competition affecting fertility leads to species 2 winning, and yet both species coexist in the model with both types of competition).

This prompted the exciting idea that community models with stage structure may foster more widespread coexistence than classic Lotka-Volterra theory, suggesting stage structure as a key missing ingredient in our explanations of diversity maintenance.

However, the results of [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF] had two limitations that prevented the generalization of this coexistence-enhancing role of stage structure. First, their investigations were restricted to a particular life cycle where juveniles and adults occupy separated habitats or niches (i.e., adults can only compete with adults and juveniles with juveniles), and these complex life cycles (sensu [START_REF] Wilbur | Complex life cycles[END_REF] were even suggested to be potentially responsible for the coexistence observed. Such segregation of competition faithfully represents some systems with metamorphosis, such as odonates (dragonflies and damselflies) or anurans (frogs), where there is an ontogenetic niche shift [START_REF] Werner | The ontogenetic niche and species interactions in size-structured populations[END_REF]. But such a model does not correspond to say, birds and mammals where adults can impose strong competition on the juveniles, as well as numerous other taxa of animals or plants for which the niche shift after development is only partial, which implies that the much larger adults can exert the bulk of the competition felt by juveniles. Thus whether competition acting on multiple life-history parameters could foster coexistence is still an open question for general life cycles, and whether complex life cycles have anything special that promotes coexistence is also something to uncover. Second, the analyses of [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF] relied entirely on numerical simulations, which means that while we can observe emergent coexistence in the Ricker-based density-dependent models that they considered, it is still difficult to understand fully the phenomenon from a theoretical standpoint. Moreover, it is yet unclear if the emergent coexistence equilibria highlighted by [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF] are structurally stable (robust to small deviations of the model framework, such as a slightly different shape of densitydependence) and realistic (occurring for a large range of parameter values). This highlights the need to investigate further structured population models with additional model structures and analytical techniques.

A further key question is whether emergent coexistence can occur in a S-species context (where S 2).

Indeed, the emergent coexistence observed by [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF] resulted from a trade-off across the life-history stages which allowed each species to the best competitor for a given vital rate. While it is obvious that such a trade-off may occur with two species, there are multiple ways to generalize such tradeoffs in competitive rankings to the S-species case, and not all of them necessarily lead to more coexistence.

In fact, unstructured Lotka-Volterra models suggest great caution in generalizing results from two-species competition to many-species competition, as the criteria for coexistence with many species usually correspond to much broader niche separation between species than required for two species, with intraspecific competition dominating [START_REF] Barabás | The effect of intra-and interspecific competition on coexistence in multispecies communities[END_REF]. Similar phenomena might be expected for stage-structured competition models.

Here, we extend the results of [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF] to a new model framework using the canonical Beverton-Holt density-dependence (which is closer to the continuous-time logistic model, [START_REF] Cushing | Some discrete competition models and the competitive exclusion principle[END_REF] and life cycles where juveniles can experience competition from adults. We present an analytical invasion criterion applicable in this new model as well as to the metamorphosis model structure numerically explored by [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF]. This criterion allows us to better understand how and why emergent coexistence occurs in both models, and is valid whenever the single-species models have stable fixed points, which occurs for iteroparous life cycles. We also use the sensitivity analysis developed by Barabás et al. (2014a) to explore the robustness of emergent coexistence equilibria. We find that emergent coexistence is a general feature of structured two-species competition models with two density-dependent vital rates, and that this equilibrium is not more sensitive to perturbation of parameters than 'classic' coexistence equilibria allowed by niche separation. We finally explore S-species models' properties using simulations of diverse community structures (number of initial species, similarity of species) and for both life cycles. These manyspecies models show the limitation of the emergent coexistence mechanism to explain coexistence in highly speciose communities.

2 Models and Methods

Competition models

Our main model is a two-species difference equation with two life stages (juvenile and adult) for each species.

It is a development of models presented in [START_REF] Fujiwara | Coexistence of competing stage-structured populations[END_REF], which were restricted to the case where only one matrix parameter per species (either fertility or survival) is affected by densities. Our model extends this framework to multiple density-dependent parameters per species, which provides a useful springboard to combine two types of competition, affecting both fertility and survival rates. Our model could also be viewed as an extension to iteroparous life-histories and delayed maturation of the model of [START_REF] Cushing | Coexistence of competing juvenile-adult structured populations[END_REF], who considered two-species competition with semelparous life-histories and fast maturation. Our model resembles [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF]'s in that two vital rates are density-dependent, and in that the life cycle includes iteroparity and a delayed maturation; however, we consider like [START_REF] Fujiwara | Coexistence of competing stage-structured populations[END_REF][START_REF] Cushing | Coexistence of competing juvenile-adult structured populations[END_REF] Beverton-Holt functions to model competition. These competition models are sometimes called Leslie-Gower after [START_REF] Leslie | The properties of a stochastic model for two competing species[END_REF], and are the closest discrete-time equivalents of continuous-time Lotka-Volterra competition models [START_REF] Cushing | Some discrete competition models and the competitive exclusion principle[END_REF]. They promote fixed points equilibria, instead of cycles or chaos that can be generated by Ricker functions at high fecundities [START_REF] Neubert | Density-dependent vital rates and their population dynamic consequences[END_REF], though see [START_REF] Cushing | Coexistence of competing juvenile-adult structured populations[END_REF] for semelparous life-histories and associated two-cycles. The Beverton-Holt functions help to find analytically the equilibria whenever multiple vital rates are density-dependent, which then fosters analytical insights with invasion criteria.

In our main model (model 1), competition is generated by adults because it felt more realistic for many animal species in which adults and juveniles share the same habitat (in mammals and birds, see e.g. snowshoe hares populations, [START_REF] Boutin | The effect of conspecifics on juvenile survival and recruitment of snowshoe hares[END_REF] or avian scavengers populations, [START_REF] Wallace | Competitive interactions within and between species in a guild of avian scavengers[END_REF].

Without any competition, population densities of the two stages of a species after one time step can be expressed (we arbitrarily pick species 1 for illustrative purposes):

n 1j (t + 1) = n 1j (t)(1 -γ 1 )s 1j + n 1a (t)f 1 (1)
and

n 1a (t + 1) = n 1j (t)γ 1 s 1j + n 1a (t)s 1a (2)
where γ 1 , f 1 , s 1j and s 1a respectively denote maturation rate, fertility, juvenile survival rate, and adult survival rate of species 1, all in [0, 1] except

f 1 ∈ R + . Additionally requiring γ i > 0 and s i,a < 1 (i ∈ {1, 2})
is necessary to ensure that individuals do not get infinite youth or infinite lifespans, respectively.

Then, we add competition on juvenile survival rates and fertility by making these parameters densitydependent. In our first and main model, we consider that competition is generated by adults of both species.

Using the Beverton-Holt model of competition, we get:

f 1 = π 1 1 + α 11 n 1a + α 12 n 2a (3)
and

s 1j = φ 1 1 + β 11 n 1a + β 12 n 2a (4)
where competition on fertility is composed of intra-specific competition (α 11 n 1a ) and inter-specific competition (α 12 n 2a ), and similarly for juvenile survival competition. Here we denoted φ 1 and π 1 for maximal juvenile survival and fertility of species 1.

With the vector n i = (n ij , n ia ), we have in matrix form

n i (t + 1) = A i (n(t))n i (t). Our model (model 1)
can be written for one species with a projection matrix A i , for instance with i = 1,

A 1 (n) =    (1-γ1)φ1 1+β11n1a+β12n2a π1 1+α11n1a+α12n2a γ1φ1 1+β11n1a+β12n2a s 1a    =    (1 -γ 1 )s 1j (n) f 1 (n) γ 1 s 1j (n) s 1a    (5)
where α ij and β ij are the competition coefficients associated to the effect of species j on species i on fertility and juvenile survival, respectively.

Although the mathematical derivations of equilibria do not impose that s i,a > 0 (iteroparity), we additionally impose that condition because the invasion and sensitivity analyses performed later require fixed point equilibria, while semelparous life-histories also produce more varied attractors, in particular 2-cycles [START_REF] Cushing | Coexistence of competing juvenile-adult structured populations[END_REF][START_REF] Cushing | Stable bifurcations in semelparous leslie models[END_REF].

We also considered a metamorphosis structure model (model 2) with a life cycle similar to the one described by [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF] and [START_REF] Cushing | Coexistence of competing juvenile-adult structured populations[END_REF], where juveniles compete only with juveniles and adults only with adults. Unlike [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF] but similarly to [START_REF] Cushing | Coexistence of competing juvenile-adult structured populations[END_REF], we chose Beverton-Holt density-dependence to model competition so that it can still be compared to model 1.

However, our model 2 also differ from [START_REF] Cushing | Coexistence of competing juvenile-adult structured populations[END_REF], who focus instead on a different life history context, with semelparity (s a = 0), fast maturation (γ = 1), and solely intraspecific competition on fertility (α 21 = α 12 = 0). In our model 2, juvenile survival does not depend on adult densities (as life-stages are separated) but depends on juvenile densities. We have a density-dependent juvenile survival rate of species 1 given by:

s M 1j (t) = φ 1 1 + β 11 n 1j (t) + β 12 n 2j (t) . ( 6
)
This model 2 can be written for one species with a projection matrix A M i , for instance with i = 1 and the same notations as for model 1:

A M 1 (n) =    (1-γ1)φ1 1+β11n1j +β12n2j π1 1+α11n1a+α12n2a γ1φ1 1+β11n1j +β12n2j s 1a    . (7)
We extended our models to S-species contexts, allowing each species to potentially interact with every other.

For model 1, we obtained the following projection matrix (for species i):

A i (n) =     (1-γi)φi 1+ S k=1 β ik n k,a πi 1+ S k=1 α ik n k,a γiφi 1+ S k=1 β ik n k,a s i,a     . (8)
The extension for model 2 is similar:

A M i (n) =     (1-γi)φi 1+ S k=1 β ik n k,j πi 1+ S k=1 α ik n k,a γiφi 1+ S k=1 β ik n k,j s i,a     . (9)

Scenarios of coexistence and associated parameter sets

We chose three scenarios of coexistence for illustration and sensitivity analysis, matching those described by [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF]: (1) a classic coexistence state that is indicated by both α and β coefficients (equivalent to the usual intra > interspecific competition coefficients in the factorised formulation of the unstructured Lotka-Volterra model, e.g. [START_REF] Letten | Linking modern coexistence theory and contemporary niche theory[END_REF], (2) an emergent coexistence state where the competition coefficients associated with the two vital rates indicate exclusion by the opposite species, and (3) a coexistence state where there is a priority effect indicated by competition coefficients acting on one vital rate while competition acting on the other vital rate indicates coexistence. For simplicity, we refer to these two last scenarios as 'emergent coexistence' since coexistence is not guaranteed by the conditions verified separately by each set of coefficients. We note that coexistence scenario (3) does not strictly correspond to an emergent outcome as defined by [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF] because competition coefficients for one of the two vital rates in fact already indicate coexistence.

Although we consider below arbitrary parameters in Invasion analysis, we also selected for numerical applications some parameter sets. To obtain parameter sets promoting emergent coexistence, we set a tradeoff between α ij and β ij for both species (large α ij is associated to low β ij and vice versa).

To be able to compare the three scenarios, we set intraspecific competition coefficients to a constant and only changed interspecific competition coefficients. Note that we also chose the values of interspecific competition coefficients such that we obtain the same densities at equilibrium across the three coexistence scenarios. The equilibrium densities have been obtained through simulation of the dynamics over 3000 time steps. No closed-form solutions for the equilibrium densities are known for all possible values of parameters, but some analytical progress can be made in special cases (Supplement A). Parameter values for the three corresponding parameter sets are given in Table 1.

Set π 1 π 2 γ 1 γ 2 φ 1 φ 2 s
= β ii = 0.1.
The three parameter sets led to the same densities at equilibrium with n 1 = 200 ± 1 and n 2 = 130 ± 1 (these have been fine-tuned by changing interspecific interactions while keeping other parameters constant).

To check what coexistence outcomes were indicated by the competition coefficients associated to fertility and juvenile survival, we used the invasion criteria given by [START_REF] Fujiwara | Coexistence of competing stage-structured populations[END_REF], valid for models with a single density-dependent vital rate affected by competition (either fertility or juvenile survival). For example, the invasion criteria of the model with only fertility affected by competition (i.e., with only α coefficients) is

given by:

R α (i) = R i R j α jj α ij ( 10
)
where

R i = π i f (c) i -1 (11)
and

f (c) i = 1 γ i s j (1 -s a ) (1 -s j + γ i s j ) . ( 12

Invasion analysis of exclusion equilibria

We devised analytical invasion criteria for our model combining competition on fertility and juvenile survival, starting with model 1 where all competition is generated by adults. These differ from the invasion criteria of the previous section in that both α and β are now potentially non-zero. We were inspired by the work of [START_REF] Cushing | An introduction to structured population dynamics[END_REF]Cushing ( , 2008) ) who studied the stability of the exclusion equilibria (with the focal invading species absent) to describe the ability of the focal species to invade a community. We therefore looked for conditions equivalent to an unstable exclusion equilibrium, which corresponds to the fact that the absent species can actually invade. All equations are presented in Appendix 1.

We used the aggregated parameters C = (1-γ)φ and D = πγφ 1-sa in the single-species model. By considering the model with a single species (since at the time of invasion, the invader is almost absent), we found an expression for the inherent net reproductive number, C + D = (1 -γ)φ + πγφ 1-sa . This number must exceed 1 to have a viable population of the resident species. We then obtained the densities at a stable positive equilibrium for the single-species model:

     n * j = (1 -s a ) 1+βn * a γφ n * a n * a = (αC-α-β)+ √ (-αC+α+β) 2 -4αβ(-C-D+1) 2αβ (13)
provided that

C + D = (1 -γ)φ + πγφ 1 -s a > 1. ( 14
)
We then evaluated the stability of the exclusion equilibrium in the two-species model: if it is locally asymptotically stable, the invader converges to zero density; if not, the invader invades. To do this, we calculated the eigenvalues of the Jacobian of the 2 species system evaluated at the exclusion equilibrium. We found the following condition for stability when species 2 is absent and species 1 is resident:

C 2 1 + β 21 n * 1a + D 2 (1 + β 21 n * 1a )(1 + α 21 n * 1a ) < 1 (15)
with C 2 = (1-γ 2 )φ 2 and D 2 = π2γ2φ2 1-s2a as denoted previously. This last expression (eq. 15) gives us an invasion criterion, in the sense that if the expression is larger than 1, species 2 is able to invade the community when rare.

We applied the exact same method to model 2 with metamorphosis structure, where juveniles compete with juveniles and adults compete with adults (Supplement B).

We found analytical expressions for the boundaries of coexistence domains over α 12 and α 21 . These expressions are respectively given by equations 16 and 17 that are derived directly from the invasion criteria:

α * 21 = D 2 (1 + β 21 n * 1a ) -C 2 -1 1 n * 1a (16) α * 12 = D 1 (1 + β 12 n * 2a ) -C 1 -1 1 n * 2a (17)
If both conditions α 12 <α * 12 and α 21 <α * 21 are met, then there is coexistence. If both conditions are violated, there is a priority effect, and if one condition is met and the other is violated, there is exclusion of one of the species.

To visualize a case of emergent coexistence outcomes as did [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF] for their Ricker competition model, we represent the coexistence domains over α 12 and α 21 with fixed values of β 12 and β 21 , for model 1. We stress that here, these coexistence domains arise from analytical expressions of the boundaries rather than numerical simulations of the dynamical systems. We computed α * 12 and α * 21 for the three parameter sets mentioned in section 2.2 where β coefficients indicated respectively coexistence, exclusion of species 1 and a priority effect.

Sensitivity analysis

To examine the sensitivity of equilibrium densities to changes in parameters under different scenarios of coexistence, we applied the community-level sensitivity analysis developed by Barabás et al. (2014a) to model 1. The computation of fixed point densities, required to perform the analysis, was done numerically.

With iteroparity, only stable fixed points were observed, both in single-species and two-species competition models, as shown in Supplement C. The two-species competition models fixed points do not have closed-form expressions in the general case-although this may sound surprising (but see [START_REF] Meszéna | Adaptive dynamics in a 2-patch environment: a toy model for allopatric and parapatric speciation[END_REF], it was indeed easier to compute invasion criteria than to find analytically the two-species fixed point, that may be analytically tractable only in special cases (see Supplement A). We sum up here the notations and expressions used to compute the sensitivity. Explaining all the steps would require paraphrasing too much of Barabás et al. (2014a), so we simply refer the reader to their paper for more details. We use n i = n i q i where n i is the vector of population density of the species i at equilibrium for the various stages and q i a vector of weights given to stages. We considered only the case where q i is a vector of ones. We also need to define p i = ni ni the proportion of each stage of species i. We denote A i the projection matrix of species i with its eigenvalues λ k i and its left and right eigenvectors v k i and w k i . We also define the dominant eigenvalue associated with its eigenvector : λ i and w i . We normalized the vectors to have |w i | = 1 and v j i w k i = δ jk with δ jk = 1 if j = k and 0 otherwise. Barabás et al. (2014a) used regulating variables R µ which correspond for example for species 1 in our model 1 to R α,1 = α 11 n 1a + α 12 n 2a and R β,1 = β 11 n 1a + β 12 n 2a . Note that the notation for regulating variables of Barabás et al. (2014a) is similar to that used for invasion criteria in single-species models that we previously introduced, even though ecological meanings differ greatly. As regulating variables are not used later in the text, there is no risk of confusion and we simply follow the conventions of Barabás et al. (2014a) in this section.

We then computed the sensitivity of population densities of both species to perturbations of each parameter using the following expressions, where E is the perturbed parameter:

dn i dE = S j=1 (M ij ) -1 g j ( 18
)
where the matrix M ij is given by

M ij = µ,ν v i ∂A i ∂R µ w i δ µν - ∂G µ ∂R ν -1 ∂R ν ∂n j p j ( 19
)
and the vector g j is given by

g j = v j ∂A j ∂E w j + σ,ρ v j ∂A j ∂R σ w j δ σρ - ∂G σ ∂R ρ -1 ∂G ρ ∂E . ( 20
)
The function G µ (R ν , E), introduced by Barabás et al. (2014a) to simplify notations, is given by

G µ (R ν , E) = S j=1 n j q j w j ∂R µ ∂n j sj k=2 1 λ j -λ k j × w k j - q j w k j q j w j w j ⊗ v k j A j (R ν , E)w j . ( 21
)
With these expressions, we were able to compute the sensitivity to perturbations of parameter values for the three parameter sets / contrasted scenarios of coexistence. These have near-equal densities at equilibrium (see section 2.2), to allow for meaningful comparisons.

S-species simulations

We simulated the dynamics of larger communities (5, 10 and 40 species) with parameters chosen to favour emergence of coexistence through life-history complexity, using the extension of model 1 and of model 2 in S-species contexts (eqs. 8 and 9). Our idea to extend the mechanism behind emergent coexistence from two to an arbitrary number of species was to create opposite competitive hierarchies in the competition coefficients related to fertility (α) and juvenile survival (β). We therefore devised a mechanism by which a species strongly competitive on fertility would be weakly competitive on juvenile survival and reciprocally.

This was done through a negative correlation between the coefficients α ij and β ij . The negative correlation was parameterized using a bivariate normal distribution and a covariance between the two variables equal to -0.9σ 2 , where σ 2 is the variance of both α ij and β ij .

To draw the full set of parameter values, we have used:

• Log-normal fertilities π i ∼ Log-N (µ π , σ 2 π )
• Beta-distributed survival and transition probabilities

φ i ∼ Beta(a φ , b φ ), s a,i ∼ Beta(a sa , b sa ) and γ i ∼ Beta(a γ , b γ ) • Normally distributed intraspecific competition coefficients α ii ∼ N (µ α , σ 2 α ) and β ii ∼ N (µ β , σ 2 β )
• As mentioned above, when i = j we use interspecific coefficients drawn as

(α ij , β ij ) ∼ N 2 (µ αβ , Σ αβ ) with Σ αβ =    σ 2 ρσ 2 ρσ 2 σ 2    and ρ = -0.9. We have µ αβ = 0.8    µ α µ β    to ensure a strong interspecific
competition, which seemed necessary to see possible effects of trade-off between competitive coefficients.

Meta-parameter values determining the abovementioned parameter distributions are given in full in Appendix 3. We chose 3 sets for the distribution-level parameters (µ π , σ π , a φ , b φ , ...), corresponding to 3 different levels of variance in the vital rates that allowed us to have either species with very close vital rates (i.e., low variance of the parameter distributions across species) or more variable vital rates (medium and high variation, see

Appendix 3).

We simulated the dynamics and computed the number of extant species after 3000 time steps, which allowed to reach an equilibrium state (see Appendix 3). Then, for each of the three main parameter sets, we created new, replicated parameter sets where we permuted all the inter-specific coefficients in order to remove the correlation between α ij and β ij , while keeping the same inter-specific competition values. We performed these permutations 200 times and represented the histogram of the number of species that coexist after 3000 time steps for the 200 permutated datasets. This allowed us to compare scenarios where α ij and β ij are negatively correlated to the null hypothesis of no correlation between them (which renders the emergence of coexistence through stage structure impossible).

Results

Invasion analysis

First of all, the analytical expressions for the boundaries of coexistence domains obtained thanks to the invasion analysis (eq. 16 and 17) can be interpreted to understand how the trade-off between competition coefficients actually promotes coexistence. In equation 16 (and similarly in equation 17), as n 1a , C 2 , D 2 are independent of inter-specific competition coefficients, a decrease of β 21 (the effect of species 1 on species 2's juvenile survival) on the right hand side will mechanically increase the value of α * 21 on the left-hand site (the boundary value for α 21 , the effect of species 1 on species 2's fertility). What does this mean ecologically?

When the effect of species 1 on species 2's juvenile survival β 21 decreases, α 21 can increase without species 2 getting extinct. This means that as one type of competition is lowered, the other type of competition can be increased while having both species coexisting. Thus we can have coexistence with low β 21 and high α 21 , and reciprocally, using the other equation, with high β 12 and low α 12 . That is, we have coexistence through opposite competitive hierarchies, which can be obtained for the full range of competition coefficients.

To better visualize coexistence properties implied by our invasion analysis, we produced coexistence outcome domains for model 1. These are similar to [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF]'s but are computed analytically.

They are produced using the newly derived invasion criterion (eq. 15). Figures 1A and1B are used for reference (with dotted lines corresponding to the domains boundaries of figure 1A, which materializes what coexistence outcomes are predicted based on α coefficients). We found, as shown in Figs. 1C, 1D and1E, that emergent coexistence outcomes can be observed in our model. This is best seen in Fig. 1D, where the coexistence domain in our model combining α and β coefficients (in light grey) is extended far to the right of the reference dashed line, well beyond what would be indicated based on α coefficients only (Fig. 1A). Thus we have coexistence with competition on both vital rates when competition on a single vital rate (α) would instead have indicated exclusion of species 2, while β indicates exclusion of species 1. The three parameter sets highlighted in the Models and Methods section are marked by asterisks, as panels C, D, E actually include a larger span of parameter space since (α 12 , α 21 ) do vary.

In Figs. 1C and1E, we do not find strictly emergent coexistence sensu [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF], as they restricted this label to coexistence which is not indicated by either α or β but is obtained by combining both types of competition. We do find other emergent outcomes though (e.g. in Fig. 1E: exclusion of one of the species when α or β would indicate priority effect and coexistence, respectively). We also reproduced Fig. 1 for model 2 where habitats of juveniles and adults are separated so that they cannot compete with each other, as in the metamorphosis structure model studied by [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF] (see Supplement B for details). The results are very similar.

Sensitivity analysis

Sensitivity and elasticity computation at the equilibrium for the three scenarios provides information relative to how much a small perturbation can change the equilibrium densities, and therefore how robust coexistence is to change in parameter values: if a small change in a parameter can make the density of a species greatly decrease, it is not robust. With some approximation (a linear one), it also provides information on how a large perturbation in parameter values could generate exclusion of one species by moving the equilibrium densities of this species to 0 (see for this the coexistence regions in Appendix 2). We represented elasticity E ni dni dE , i.e., a relative response of the density n i of species i to a relative perturbation on E. The elasticity values obtained from the sensitivity analysis performed on our three parameter sets are given in Fig. 2. 2.

-2 0 2 p 1 p 2 g 1 g 2 f 1 f 2 s 1a s
We found here that the emergent coexistence that occurs when there are opposite competitive hierarchies on fertility and juvenile survival (parameter set 2) is as sensitive as the 'classical' (niche-based) coexistence that occurs when the competition coefficients on both vital rates indicate coexistence (parameter set 1).

However, when coexistence stems from competition on one vital rate indicating a priority effect and the other vital rate indicating niche-based coexistence (parameter set 3), the resulting equilibrium is more sensitive to perturbation, according to our elasticity analysis. Note that sensitivity/elasticity values are formally valid only close to the fixed point equilibrium. The coexistence regions representation (Barabás et al. 2014b, presented in Appendix 2), which stretches a bit that validity, confirmed the results of Fig. 2 by showing similar regions for the 'classical' and the 'emergent' coexistence, together with a smaller domain for the third scenario (priority effect on β + coexistence on α → coexistence). This is likely due to the destabilizing influence of priority effects.

S-species simulations

The species richnesses at the end of the simulations for the extended model 1 with S species (see eq 8) are

given in Fig. 3 (see Supplement D for results with extended model 2). Each panel displays a histogram of the number of species remaining in the community after 3000 time steps, for the original parameter set including a negative correlation between α ij and β ij , as well as for the 200 permutated parameter sets where this correlation is set to zero.

species, low variance

Number of species at equilibrium We checked that the negative correlation between α ij and β ij indeed induces opposite hierarchies of competition in sets of coefficients associated to fertility and juvenile survival (i.e., a negative correlation between rank(α ij ) and rank(β ij ), see Appendix 3). In spite of this, the parameter sets with the above-mentioned opposite competitive hierarchies do not systematically lead to richer communities than the corresponding permutated parameter sets (where Corr(α ij , β ij ) = 0). For very close parameters between species (low variances in parameter distributions, Figs. 3A, 3B and 3C), opposite competitive hierarchies seem to have a positive effect on the species richness, when compared to the null hypothesis of zero correlation between α ij and β ij . However, the assumption of a community composed solely of very similar species for all life-history parameters (as opposed to similar averaged birth and death rates, as in neutral theory) is rather unrealistic.

To sum up, we showed that the negative correlation between α ij and β ij , even if it indeed induces opposite competitive hierarchies on the two vital rates, does not seem to allow, for reasonably variable distributions of parameter values, an increase of the number of species coexisting at the equilibrium. In other words, two-species coexistence through differing competitive hierarchies on life-history parameters does not scale up to many-species.

Discussion

In this article, we studied competitive interactions between species with structured life-cycles, with multiple vital rates that can depend on population densities. Using structured models for two competing species, [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF] described a form of coexistence whereby opposite competitive hierarchies in the competition coefficients associated to two vital rates could yield coexistence in the model where these two density-dependent vital rates combine-while models with a single density-dependent vital rate would predict exclusion of the inferior competitor. To show the existence of such emergent coexistence, they used Ricker functions for density-dependence and a particular interaction setup where competition between adults and juveniles is forbidden, as in some systems with metamorphosis, where adults and juveniles occupy separated habitats. This brought into question whether these results could be generalized to other functional forms and life cycles. We considered as a baseline an alternative and perhaps more frequent life cycle setup, where adults can affect both adults and juveniles, so that both fertilities and juvenile survival rates can depend upon adult densities of both species. We first found, using structured two-species models, that emergent coexistence can be generalized to Beverton-Holt functional forms and our new life cycle, in addition to the life cycle with metamorphosis considered by [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF] (see Supplement B). We provided an in-depth explanation of how emergent coexistence can occur, using analytical invasion criteria (facilitated by the use of Beverton-Holt functions). Then, we explored the robustness of such emergent coexistence equilibria with a community-level sensitivity analysis and we found that a perturbation on parameters had not more effect on equilibrium densities in the case of emergent coexistence than in a 'classical' coexistence allowed by niche separation. Finally, we numerically explored the possibility for such mechanisms of emergent coexistence to manifest in S-species models (with S >> 2). We observed that opposite competitive hierarchies acting on fertility vs juvenile survival, allowed by a negative correlation between the two sets of competition coefficients, did not seem to generate substantial increases in species richness in our simulated communities. Importantly, this was true both for our model 1 where all competition is generated by adults, as well as model 2 where adults only compete with adults and juveniles with juveniles-had [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF] considered more species they would have found similar results to us.

We now dive deeper into each of the abovementioned results. As we succeeded to reproduce the emergent coexistence of [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF] with a general invasion framework and an additional, frequently observed life cycle, we think that coexistence emerging through opposite competitive hierarchies on different life history traits may be a general feature of two-species two-life-stage models. We found nothing in complex life histories with habitat/niche shifts that specifically promotes this effect of stage structure on coexistence.

This reinforces the importance of structured models to study interactions between species and their role in shaping community dynamics [START_REF] Miller | Thinking inside the box: community-level consequences of stagestructured populations[END_REF], as we find that the set of life cycles where stage structure affects coexistence is considerably extended. Furthermore, by finding an analytic invasion criterion, we ensured that emergent coexistence is not dependent on specific sets of vital rate parameters but rather a general feature of two-species two-stages models, depending on competition hierarchies. The robustness of such coexistence equilibria to perturbation of parameters has then been confirmed by the communitylevel sensitivity analysis, as the emergent coexistence scenario was not more sensitive than the classical coexistence scenario (where coexistence was indicated by α and β coefficients separately). To our knowledge, there has been little implementation in practice of the discrete-time sensitivity analysis for coupled nonlinear systems developed by Barabás et al. (2014a). Here, it has allowed us to explore in more depth the properties of equilibria, as we could not easily find an expression for equilibrium densities (they are defined by the intersection of two conic sections as shown in Supplement A, and parameters are so interwoven in those expressions that even if analytical solutions were to be found, these might be too complex to interpret). The robustness of the emergent coexistence equilibria to perturbations further suggests that such coexistence may be likely even when environmental stochasticity is added to two-species two-stages models. It is therefore an explanation to keep in mind when seeing coexisting pairs of species in the wild in spite of coexistence hierarchies on a single resource axis or vital rate indicating exclusion.

However, two-species coexistence, even with many stages, can be quite different from many-species coexistence. Our study of S-species communities did not show a positive effect of opposite competitive hierarchies in the two sets of competition coefficients, except in communities where species were very similar in all their parameters-but these communities are probably not very lifelike. We have done additional explorations of many-species models, using slightly different correlation structures, such as negative correlations between interspecific α and β coefficients within each pair of species (see Appendix 3), rather than simply a negative correlation between pooled α ij and β ij values for all pairs of species. This stronger condition on negative correlations between competitiveness on the two vital rates had no effect on the species richness of the fi-nal community. In sum, we showed that the properties promoting emergence of coexistence through stage structure in the 2-species models could hardly be extended to many-species models where all species interact together (albeit sometimes weakly). These results resonate with those of simpler, unstructured models [START_REF] Barabás | The effect of intra-and interspecific competition on coexistence in multispecies communities[END_REF]: criteria needed to promote S-species coexistence with S 2 are much more stringent than those for two species, and require a more important separation of niches (manifesting in the ratio between intra and interspecific competition coefficients).

Finally, we should keep in mind that even if our results are both robust and supported by analytical formulas, some assumptions could be debated. While we considered iteroparous life histories, Cushing et al.

(2007) considered semelparous life histories with more variable dynamics, which creates additional possibilities for coexistence, with 2-cycle attractors where juveniles of species 1 co-occur with adults of species 2 and vice versa. The extension of these semelparous life histories to our new life cycle would be quite interesting. More generally, [START_REF] Cushing | Coexistence of competing juvenile-adult structured populations[END_REF], [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF] and ourselves considered a modified Lotka-Volterra framework, with Lotka-Volterra competition coefficients modulating the effect of adult or juvenile densities on several vital rates. However, we have no mechanistic derivations for such models (a derivation through MacArthur-style consumer-resource model with separation of time scales could perhaps be attempted). Structured mechanistic models of resource-based competition or competition for space might inform on the realism of the α ij and β ij correlation structures considered here [START_REF] Goldberg | Competitive effect and response: hierarchies and correlated traits in the early stages of competition[END_REF][START_REF] Loreau | Competitive exclusion and coexistence of species with complex life cycles[END_REF][START_REF] Tilman | Resource Competition and Community Structure[END_REF][START_REF] Qi | Coexistence under hierarchical resource exploitation: the role of r*-preemption tradeoff[END_REF]. That said, we expect that even with a mechanistic perspective, our main finding will likely be robust: while coexistence emerging from life-history complexity is quite easily seen in two-species pairs, it is much more difficult to obtain in many-species models. In speciose communities, some degree of niche separation is required to prevent species to interact too much with each other. This could be promoted by many weak-few strong interactions, possibly combined with some interaction network structures for each stage (e.g. [START_REF] Kinlock | Uncovering structural features that underlie coexistence in an invaded woody plant community with interaction networks at multiple life stages[END_REF] or density-dependent vital rate, as the network structure is likely to be nonrandom and could differ among stages.

authors contributed equally to subsequent versions.

Appendix 1 Invasion analysis

We develop here the method to obtain invasion criteria in two-species two-life stages models including competition on two vital rates, such as model 1 combining competition on fertility and juvenile survival.

We first place ourselves at the exclusion equilibrium where species 1 lives without species 2:

ñ =          n * 1j n * 1a 0 0          . (A1)
To compute the equilibrium densities n * 1j and n * 1a , we have to solve the following system corresponding to the model with a single species (we dropped here species-specific indices for simplicity):

     n j = (1 -γ)s j (n)n j + f (n)n a n a = γs j (n)n j + s a n a (A2)
with fertility f (n) and juvenile survival s j (n) affected by intra-specific competition only, according to the relations

f (n) = π 1 + αn a and s j (n) = φ 1 + βn a . ( A3
)
We now express n j according to n a :

n j = (1 -s a )n a 1 γs j (n) (A4) = (1 -s a ) 1 + βn a γφ n a (A5)
and then, as we have by assumption n a = 0 and n j = 0, we obtain

n j = (1 -γ)s j (n)n j + f (n)n a ⇔ 1 = (1 -γ) φ 1 + βn a + π 1 + αn a n a n j ⇔ 1 = (1 -γ) φ 1 + βn a + πφγ (1 -s a )(1 + βn a )(1 + αn a ) ⇔ (1 + βn a )(1 + αn a ) = (1 -γ)φ(1 + αn a ) + πγφ 1 -s a
that leads to the polynomial of degree 2:

αβn 2 a + (α + β -α(1 -γ)φ)n a + 1 -(1 -γ)φ - πγφ 1 -s a = 0. ( A6
)
We denote C = (1 -γ)φ and D = πγφ 1-sa and we find the solutions of the polynomial equation:

n * a = (αC -α -β) ± (-αC + α + β) 2 -4αβ(-C -D + 1) 2αβ . (A7)
We search for a positive solution and we have C = (1 -γ)φ < 1 which implies αC -α -β < 0 assuming positive interaction coefficients (α > 0 and β > 0). Positivity is then equivalent to:

(αC -α -β) + (-αC + α + β) 2 -4αβ(-C -D + 1) > 0 (A8) ⇐⇒ -4αβ(-C -D + 1) > 0 (A9) ⇐⇒ C + D > 1. ( A10
)
The equivalence is provided by the fact that the left term has its square inside the square root. We define the inherent net reproductive number as C + D = (1 -γ)φ + πγφ 1-sa that must exceed 1 to have a viable species. We finally obtain the densities at a stable positive equilibrium for the model with a single species:

     n * j = (1 -s a )( 1+βn * a γφ )n * a n * a = (α(1-γ)φ-α-β)+ (-α(1-γ)φ+α+β) 2 -4αβ(-(1-γ)φ-πγφ 1-sa +1) 2αβ (A11) provided that C + D = (1 -γ)φ + πγφ 1 -s a > 1. (A12)
We have then obtained the densities at the exclusion equilibrium. We have now to evaluate the stability of the equilibrium (n * j , n * a , 0, 0): if it is locally asymptotically stable, the excluded species cannot invade the community starting from very small density. Conversely, if it is unstable, this means that the excluded species can invade provided that equilibria of the single-species model are restricted to stable fixed points. This last restriction is important. Indeed, if the single-species model has a stable two-cycle or an even more complex fluctuating attractor, the fixed point will already be unstable, and a proper invasion analysis should instead be done on this more complex attractor. Two-cycles tend to occur for s a = 0 [START_REF] Cushing | Coexistence of competing juvenile-adult structured populations[END_REF], as shown in Supplement C. So long as s a > 0 this model does not show 2-cycles-so far as we have seen numerically and following the argumentation of [START_REF] Cushing | Coexistence of competing juvenile-adult structured populations[END_REF] for model 2. The combination s a = 0 and γ = 1 renders the projection matrix primitive, which in turns creates possibilities for bifurcations towards twocycles directly from the trivial single-species (0,0) equilibrium [START_REF] Cushing | On ebenman's model for the dynamics of a population with competing juveniles and adults[END_REF][START_REF] Cushing | Nonlinear semelparous leslie models[END_REF]. We therefore restrict ourselve to the case s a > 0, although the equations derived above in this Appendix for the single-species fixed point are also valid whenever s a = 0.

To evaluate the (un)stability of the exclusion equilibria in the two-species system, we have to compute the eigenvalues of the Jacobian of the 2 species system evaluated at the exclusion equilibrium. The full system iterated over a time step reads:

                 n 1j (t + 1) = (1 -γ 1 ) φ1 1+β11n1a(t)+β12n2a(t) n 1j (t) + π1 1+α11n1a(t)+α12n2a(t) n 1a (t) n 1a (t + 1) = γ 1 φ1 1+β11n1a(t)+β12n2a(t) n 1j (t) + s 1a n 1a (t) n 2j (t + 1) = (1 -γ 2 ) φ2 1+β22n2a(t)+β21n1a(t) n 2j (t) + π2 1+α22n2a(t)+α21n1a(t) n 2a (t) n 2a (t + 1) = γ 2 φ2 1+β22n2a(t)+β21n1a(t) n 2j (t) + s 2a n 2a (t). (A13)
We place ourselves at the abovementioned exclusion equilibrium where species 1 dominates the community and species 2 is absent. The Jacobian evaluated at this point is a 4 × 4 matrix and it has the following triangular block form (see also Cushing 2008):

   B 1 B 2 0 B 3    . (A14)
Therefore, we only need to know the eigenvalues of the 2 × 2 matrices B 1 and B 3 . The B 1 matrix corresponds to the Jacobian of the system in the absence of species 2, so provided that the single-species equilibrium is stable, we can infer that it has a spectral radius lower than 1. Unfortunately, although all simulations that we have done, even for very large fertilities or low competition coefficients, provide a stable equilibrium for C + D > 1, it is not straightforward to show local stability analytically. Indeed matrix B 1 writes:

   (1 -γ) φ 1+βna 1-γ γ -βna 1+βna (1 -s a ) + π 1+αna 1 -αna 1+αna γ φ 1+βna -βna 1+βna + s a    (A15)
with n a given by eq. ( A11), making the use of Jury conditions (used below on the simpler B 3 block) difficult here, although the determinant simplifies somewhat. We therefore conjecture the stability of the single-species equilibrium, with a spectral radius of B 1 below 1, which we deem reasonable in this context.

The B 3 matrix is given by:

B 3 =    (1 -γ 2 ) φ2 1+β21n * 1a π2 1+α21n * 1a γ 2 φ2 1+β21n * 1a s 2a    . ( A16
)
The B 3 matrix block, which is part of the Jacobian, also corresponds to the projection matrix of species 2 when invading species 1, which implies that the exclusion equilibrium is not stable when species 2 has a growth rate larger than 1 (i.e., the modulus of B 3 's leading eigenvalue is larger than 1). The stability of this two-dimensional system can be investigated thanks to the Jury conditions, given by the equations (A17, A18, A19):

1 -tr J(n) + det J(n) > 0 (A17) 1 + tr J(n) + det J(n) > 0 (A18) 1 -det J(n) > 0 (A19)
with J(n) the Jacobian matrix of the two-compartment system whose stability is being evaluated. If these three conditions hold, the system is stable. If the first condition (equation (A17)) is violated, one of the eigenvalues of J(n) is larger than 1, which means in our context that the excluded species can invade the community. The violations of the conditions (A18) and (A19) correspond to the creation of limit cycles or invariant loops [START_REF] Neubert | Density-dependent vital rates and their population dynamic consequences[END_REF]. We then check these conditions on the B 3 matrix. We have:

tr B 3 = (1 -γ 2 ) φ 2 1 + β 21 n * 1a + s 2a det B 3 = (1 -γ 2 ) φ 2 1 + β 21 n * 1a s 2a -γ 2 φ 2 1 + β 21 n * 1a π 2 1 + α 21 n * 1a .
We start with the third condition (A19):

1 -(1 -γ 2 ) φ 2 1 + β 21 n * 1a s 2a + γ 2 φ 2 1 + β 21 n * 1a π 2 1 + α 21 n * 1a > 0. (A20)
This equation will always be satisfied with biologically meaningful parameters and strictly competitive interaction between species:

• The third term of (A20) is positive because each parameter is positive,

• (1 -γ 2 ) φ2 1+β21n * 1a s 2a < 1 because 0 ≤ γ ≤ 1, 0 ≤ φ ≤ 1
and all β's are positive.

Since tr B 3 > 0, the first Jury's condition implies the second ((S29) ⇒ (S30)). Then, we only have to check the first Jury condition to describe the stability of the system:

1 -(1 -γ 2 ) φ 2 1 + β 21 n * 1a + s 2a + (1 -γ 2 ) φ 2 1 + β 21 n * 1a s 2a -γ 2 φ 2 1 + β 21 n * 1a π 2 1 + α 21 n * 1a > 0 (1 -γ 2 ) φ 2 1 + β 21 n * 1a -(1 -γ 2 ) φ 2 1 + β 21 n * 1a s 2a + γ 2 φ 2 1 + β 21 n * 1a π 2 1 + α 21 n * 1a < 1 -s 2a (1 -γ 2 )φ 2 (1 -s 2a ) < (1 -s 2a )(1 + β 21 n * 1a ) - γ 2 φ 2 π 2 1 + α 21 n * 1a (1 -γ 2 )φ 2 1 1 + β 21 n * 1a + γ 2 φ 2 π 2 1 -s 2a 1 (1 + β 21 n * 1a )(1 + α 21 n * 1a ) < 1.
We finally have the following condition for the stability of the exclusion equilibrium:

C 2 1 + β 21 n * 1a + D 2 (1 + β 21 n * 1a )(1 + α 21 n * 1a ) < 1 (A21)
g 1 g 2 f 1 f 2 s 1a s 2a
Maturation rate, juvenile and adult survival

g 1 g 2 f 1 f 2 s 1a s 2a
Maturation rate, juvenile and adult survival

g 1 g 2 f 1 f 2 s 1a s 2a
Maturation rate, juvenile and adult survival We see here smaller coexistence regions in the third scenario, consistent with this equilibrium being more sensitive (less robust) to perturbations. The population dynamics highlighted in Fig. A2 demonstrates that 3000 time steps are enough to reach a stable number of coexisting species.

We verified our hypothesis that a negative correlation between α ij and β ij lead to opposite hierarchies in sets of coefficients associated to fertility and juvenile survival by plotting the rank of α ij against the rank of β ij . For 40 species, these ranks are plotted in Fig. A3. The negative correlation between ranks in Fig. A3 confirmed that the negative correlation between α ij and β ij was creating the expected opposite competitive hierarchies.

However, opposite competitive hierarchies may not be sufficient to promote reciprocal exclusion (e.g., exclusion of 1 by 2 in the fertility competition model and of 2 by 1 in the juvenile survival competition model) of species within each or most pairs of species. To make the comparison between 2 and S-species model, we have to check that we consider situations that are fully comparable, i.e., where species would exclude each other out in a two-species contest on either α or β competition. We therefore checked if the trade-off between being competitive on α and β could promote such reciprocal exclusion or priority effects (in models with a single parameter that is density dependent) within pairs of species. We computed the invasion criteria of such simple models (where only one vital rate is affected by competition) for each pair of species composing the communities. For a community of initially 40 species and a small variance in parameters across species, we found the coexistence outcomes in simpler models highlighted in Fig. A4.

We found that for a community of 40 species, the proportions of pairs of species where a priority effect or a reciprocal exclusion are indicated (by the simpler models with density-dependence on a single vital rate) are low when we assume only a negative correlation between α ij and β ij . In fact, when we simply use a negative correlation to create opposite competitive hierarchies, we draw our couple of values (α ij , β ij ) in a Gaussian distribution, and then most of the values are close to their means, which does not allow to obtain pairs of species where each species is strongly competitive on a different vital rate, and exclude the other if the competition was on this vital rate only. Therefore, even if there are opposite competitive hierarchies within pairs of species modelled, this setup is unlikely to generate situations where species exclude each other when considering single-vital-rate-competition.

To consider S-species scenarios where species can exclude each other out if competition was solely on α or solely on β (hereafter referred to as 'reciprocal exclusion'), we changed our structure of correlation to generate situations where within most pairs of species, each species is strongly competitive on a different vital rate (e.g. species 1 excludes species 2 in the model with only fertility competition and vice-versa in the model with only juvenile survival competition).

Our method to achieve reciprocal exclusion of pairs (for single-vital-rate-competition) consisted in drawing pairs (a ij , b ij ) from a bivariate normal distribution with µ aij <µ bij and Corr(a ij , b ij ) < 0. We then always assign (a ij , b ij ) to (α ij , β ij ) and (b ji , a ji ) to (α ji , β ji ). In this way, for each pair of species (i, j), species i is competitive on α and species j is competitive on β. We acknowledge that this may be hard to grasp, and suggest to the reader to use the code to create an example.

We computed the invasion criteria for each pair of species and counted how often each invasion scenario in single-vital-rate-competition models, which we represented in the histogram of Fig. A5. We succeeded to generate reciprocal exclusion within most pairs of species. Moreover, we showed that a permutation of competition coefficient was sufficient to remove the structure of reciprocal exclusion (Fig. A6). However, as for our simple negative correlation between competitive ranks for α and β, this new parameter set where reciprocal exclusion occurs within most pairs of species does not allow to significantly increase the number of extant species at t = 3000 time steps (Fig. A7). Therefore, even if competition is highly structured in a way that makes species exclude each other out when a single vital rate is density-dependent, which should greatly promote situations of emergent coexistence, we did not observe a positive effect of this structure on species richness in a community larger than S = 2. We therefore conclude that emergent coexistence is highly unlikely in many-species contexts.

species, low variance

Supplement A Fixed points for the two-species models: equations

and analytical progress in a simple case

To compute the densities of both species at equilibrium for model 1, we have to solve the following system:

                 n 1j = (1 -γ 1 ) φ1 1+β11n1a+β12n2a n 1j + π1 1+α11n1a+α12n2a n 1a n 1a = γ 1 φ1 1+β11n1a+β12n2a n 1j + s 1a n 1a n 2j = (1 -γ 2 ) φ2 1+β22n2a+β21n1a n 2j + π2 1+α22n2a+α21n1a n 2a n 2a = γ 2 φ2 1+β22n2a+β21n1a n 2j + s 2a n 2a . (S1)
We start by expressing n 1j according to n 1a :

n 1j 1 - (1 -γ 1 )φ 1 1 + β 11 n 1a + β 12 n 2a = π 1 1 + α 11 n 1a + α 12 n 2a n 1j = (1 + β 11 n 1a + β 12 n 2a )π 1 (1 + α 11 n 1a + α 12 n 2a )(1 + β 11 n 1a + β 12 n 2a -(1 -γ 1 )φ 1 n 1a .
Then, we insert the latter expression in the expression for n 1a :

n 1a = γ 1 φ 1 1 + β 11 n 1a + β 12 n 2a (1 + β 11 n 1a + β 12 n 2a )π 1 (1 + α 11 n 1a + α 12 n 2a )(1 + β 11 n 1a + β 12 n 2a -(1 -γ 1 )φ 1 n 1a + s 1a n 1a .
We place ourselves in the case where both species coexist. Therefore, n 1a = 0 and we have:

(1 + α 11 n 1a + α 12 n 2a )(1 + β 11 n 1a + β 12 n 2a -(1 -γ 1 )φ 1 ) = γ 1 φ 1 π 1 1 -s 1a (1 + α 11 n 1a + α 12 n 2a )(1 + β 11 n 1a + β 12 n 2a -C 1 ) = D 1 C 1 1 + β 11 n 1a + β 12 n 2a + D 1 (1 + β 11 n 1a + β 12 n 2a )(1 + α 11 n 1a + α 12 n 2a ) = 1.
This last expression correspond to the growth rate of the species 1 that must be equal to 1 at the equilibrium.

By symmetry, we have the same equation for the other species :

(1 + α 22 n 2a + α 21 n 1a )(1 + β 22 n 2a + β 21 n 1a -C 2 ) = D 2 C 2 1 + β 22 n 2a + β 21 n 1a + D 2 (1 + β 22 n 2a + β 21 n 1a )(1 + α 22 n 2a + α 21 n 1a ) = 1.
Finally, to compute the equilibrium densities, we have to solve the following system for n 1a and n 2a :

     (1 + α 11 n 1a + α 12 n 2a )(1 + β 11 n 1a + β 12 n 2a -C 1 ) = D 1 (1 + α 22 n 2a + α 21 n 1a )(1 + β 22 n 2a + β 21 n 1a -C 2 ) = D 2 (S2)
which can also be written in the following form with x = n 1a and y = n 2a :

    
)x + (1 -C 2 -D 2 ) = 0. (S3)
This system corresponds to two bivariate quadratics equations (quadrics). Therefore, the solutions of the system are the intersection points of the two conic sections described by each equation. These conic sections are themselves the intersection of the surface of a cone with a plane. We can classify the conic sections by looking at the discriminant of the matrix representing the conic sections.

With the standard notation of a quadratic equation:

Ax 2 + Bxy + Cy 2 + Dx + Ey + F = 0. (S4)
We use the usual matrix representation of the conic sections that have the following general form:

x T A Q x (S5) with x =       x y 1       (S6)
and

A Q =       A B/2 D/2 B/2 C E/2 D/2 E/2 F       (S7)
the matrix of the quadratic equation. Then, we denote the matrix of the quadratic form A 33 :

A 33 =    A B/2 B/2 C    (S8)
which is obtained by removing the 3rd column and the 3rd line to the A Q matrix.

With these standard notations, we know that the conic section is degenerated if det A Q = 0 and is nondegenerated otherwise. Then, if the conic section is non-degenerated, we have the following classification :

• the section is a hyperbola if and only if det A 33 < 0

• the section is a parabola if and only if det A 33 = 0

• the section is an ellipse if and only if det A 33 > 0

If the conic section is degenerated (det A Q = 0), the conic section corresponds to:

• 2 intersecting lines if and only if det A 33 < 0

• 2 parallel straight lines if and only if det A 33 = 0

• A single point if and only if det A 33 > 0 Firstly, we look at the determinant of the A Q1 matrix for the first equation of the system, corresponding to the first species:

A Q1 =       α 11 β 11 α11β12+α12β11 2 β11+α11(1-C1) 2 α11β12+α12β11 2 α 12 β 12 β12+α12(1-C1) 2 β11+α11(1-C1) 2 β12+α12(1-C1) 2 1 -C 1 -D 1       . ( S9
)
We can easily find an expression of the determinant of the matrix:

det A Q1 = (α 11 β 12 -α 12 β 11 )D 1 . (S10)
We have D 1 = 0, therefore we have det A Q1 = 0 ⇔ α 11 β 12 = α 12 β 11 .

We place ourselves in the case of α 11 β 12 = α 12 β 11 . For the first equation of the system, we have :

A 33,1 =    α 11 β 11 α11β12+α12β11 2 α11β12+α12β11 2 α 12 β 12    . (S11)
We can see that the determinant of the matrix A 33,1 is equal to 0 provided α 11 β 12 = α 12 β 11 . Then we find that the conic section corresponds to 2 parallel straight lines (a degenerate parabola) when α 11 β 12 = α 12 β 11 .

In the case where α 11 β 12 = α 12 β 11 , the conic section is non-degenerated, we have to look at the sign of det A 33 . For the first equation, we have:

det A 33,1 = α 11 β 11 α 12 β 12 -( α 11 β 12 + α 12 β 11 2 ) 2 = 1 4 (α 11 β 12 -α 12 β 11 ) 2 .
The discriminant is then always positive, which indicates that the conic sections are hyperbola. As the system is symmetric, we have the same conclusions for species 2.

The analytical computation of the intersection of the conic sections, in the general case, requires a tedious analysis that consists in degenerating one conic section to find more easily the expression of the intersection.

Notably, it requires solving of a third degree equation that would lead to analytical results which are probably so complex that they might be impossible to interpret. We have not found a way to prove the uniqueness of a strictly positive intersection point but our graphical observations on multiple parameter sets seemed to show that we can find only one or zero positive point of intersection.

In order to go further on the analysis of the intersection point, we placed ourselves in the case where α = β. This case corresponds to a model where adult and juvenile parameters react in the same way to increases in adult densities. Even if it can be biologically inaccurate, it provides us a simple framework to understand how the conic sections can intersect. With α = β, we fulfill the condition to have degenerated conic sections, therefore the computation of the intersection points is easier. The equation of the degenerated conic section for species 1 becomes in this case:

α 2 11 x 2 + 2α 11 α 12 xy + α 2 12 y 2 + α 11 (2 -C 1 )x + α 12 (2 -C 1 )y + (1 -C 1 -D 1 ) = 0. (S12)
By separating the terms in x and y, we can find the following expressions :

y = - α 11 α 12 x + α 11 α 12 1 4 C 2 1 + D 1 -1 + C 1 2 (S13)
and

y = - α 11 α 12 x + α 11 α 12 - 1 4 C 2 1 + D 1 -1 + C 1 2 (S14)
which correspond to two parallel lines. We have symmetrical equations for species 2, which allow us to compute the intersections between the two pairs of lines. We find that for both species, only one of the two lines cross the section where x > 0 and y > 0. Therefore, there is at maximum one intersection point which has both coordinates positive. This point has for coordinates:

x = α 12 α 21 α 22 α 12 -α 11 α 21 α 22 α 21 1 4 C 2 2 + D 2 -1 + C 2 2 - α 11 α 12 1 4 C 2 1 + D 1 -1 + C 1 2 (S15) y = α 11 α 22 α 22 α 12 -α 11 α 21 1 4 C 2 1 + D 1 + C 1 2 - 1 4 C 2 2 + D 2 - C 2 2 . ( S16
)
We can search the conditions to have x > 0 and y > 0.

We assume α 22 α 12 > α 11 α 21 .

The condition to have y > 0 is:

1 4 C 2 1 + D 1 + C 1 2 > 1 4 C 2 2 + D 2 + C 2 2 (S17)
and to have x > 0, we need:

α 22 α 21 1 4 C 2 2 + D 2 -1 + C 2 2 > α 11 α 12 1 4 C 2 1 + D 1 -1 + C 1 2 (S18) α 22 α 12 1 4 C 2 2 + D 2 -1 + C 2 2 > α 11 α 21 1 4 C 2 1 + D 1 -1 + C 1 2 (S19)
which is valid provided the condition to have y > 0.

Therefore, we obtained the conditions to have a fixed point whose coordinates are both positive, which correspond to the conditions of coexistence, for the model where α = β.

The search for fixed points in our model combining density-dependence on fertility and juvenile survival have allowed us to better understand how a strictly positive equilibrium (interior solution) can occur. With the simpler model where α = β, we found additional conditions for coexistence than with the invasion analysis.

However, the analytical results remain difficult to interpret and/or restricted to the simplest models.

Supplement B Invasion analysis for a model with metamorphosis life history structure

We performed our invasion analysis on the model with a life-history structure similar to [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF] and [START_REF] Cushing | Coexistence of competing juvenile-adult structured populations[END_REF], where juveniles and adults are separated, i.e. juveniles are in competition with juveniles and adults with adults. The projection matrix of this model is given in equation (S36). The exact same method as in section 2.3 was used to compute the invasion criteria with this life-history structure.

In the metamorphosis model, with our notations and Beverton-Holt competition functions, the model can be written, for species 1:

A M 1 =    (1 -γ 1 ) φ1 1+β11n1j +β12n2j π1 1+α11n1a+α12n2a γ 1 φ1 1+β11n1j +β12n2j s 1a    . (S20)
We applied the invasion analysis as for our combined model. Firstly, we search the densities of exclusion equilibria. Then, we have the following system for the model with a single species (we drop here the time indices for clarity and because the model is evaluated at the fixed point):

     n j = (1 -γ) φ 1+βnj n j + π 1+αna n a n a = γ φ 1+βnj n j + s a n a . (S21)
We can easily express n j according to n a :

n a = γφ 1 + βn j n j + s a n a ⇔ γφn j = n a (1 -s a ) + n a (1 -s a )βn j ⇔ (γφ -n a (1 -s a )β)n j = n a (1 -s a ) so that n j = (1 -s a )n a γφ -(1 -s a )βn a . ( S22
)
We replace n j by its expression in the first equation of the system (S21), that leads to the polynomial of degree 2:

Eαβn 2 a + (α(1 -C) + β(π + E)) n a + 1 -C -D (S23)
where

C = φ(1 -γ), D = πγφ 1-sa and E = (1-γ)(1-sa) γ .
We find the same condition for inherent net reproductive number than for our model 1. It is given by C + D > 1 that must be true to have a viable species. This is due to the similarity of the equations to those of model 1, which mainly differ in that factor E now appears in front of the degree 2 term. We then have positive densities for the exclusion equilibrium where one species is absent:

     n * j = (1-sa)n * a γφ-(1-sa)βn * a n * a = -(α(1-C)+β(π+E))+ √ (α(1-C)+β(π+E)) 2 -4Eαβ(1-C-D) 2Eαβ (S24) provided C + D = (1 -γ)φ + πγφ 1 -s a > 1. (S25)
We have then obtained the densities at the exclusion equilibrium. We have now to evaluate the stability of this equilibrium: if it is locally asymptotically stable, the excluded species cannot invade the community starting from very small density. Conversely, if it is unstable, this means that the excluded species can invade provided that equilibria of the single-species model are restricted to stable fixed points. This last restriction is key. Indeed, if the single-species model has a stable two-cycle or an even more complex fluctuating attractor, the fixed point will already be unstable, and a proper invasion analysis should instead be done on this more complex attractor. Two-cycles tend to occur for s a = 0 [START_REF] Cushing | Coexistence of competing juvenile-adult structured populations[END_REF], as shown in Supplement C. So long as s a > 0 this model does not show 2-cycles-so far as we have seen numerically and following the argumentation of [START_REF] Cushing | Coexistence of competing juvenile-adult structured populations[END_REF] for model 2. We therefore restrict ourselves to the case s a > 0, although the equations derived above in this Appendix for the single-species fixed point are also valid whenever s a = 0.

To evaluate the (un)stability of the exclusion equilibria in the two-species system, we have to compute the eigenvalues of the Jacobian of the 2 species system evaluated at the exclusion equilibrium. The full system iterated over a time step reads:

                 n 1j (t + 1) = (1 -γ 1 ) φ1 1+β11n1j (t)+β12n2j (t) n 1j (t) + π1 1+α11n1a(t)+α12n2a(t) n 1a (t) n 1a (t + 1) = γ 1 φ1 1+β11n1j (t)+β12n2j (t) n 1j (t) + s 1a n 1a (t) n 2j (t + 1) = (1 -γ 2 ) φ2 1+β22n2j (t)+β21n1j (t) n 2j (t) + π2 1+α22n2a(t)+α21n1a(t) n 2a (t) n 2a (t + 1) = γ 2 φ2 1+β22n2j (t)+β21n1j (t) n 2j (t) + s 2a n 2a (t). (S26)
We place ourselves at the abovementioned exclusion equilibrium where species 1 dominates the community and species 2 is absent. The Jacobian evaluated at this point is a 4 × 4 matrix and it has the following triangular block form:

   B 1 B 2 0 B 3    . (S27)
Therefore, we only need to know the eigenvalues of the 2 × 2 matrices B 1 and B 3 . The B 1 matrix corresponds to the Jacobian of the single-species dynamical system, which has always been found to be stable numerically, thus for reasons explained in more detail in Appendix 1 we assume the eigenvalues have a modulus smaller than 1. The B 3 matrix is given by:

B 3 =    (1 -γ 2 ) φ2 1+β21n * 1j π2 1+α21n * 1a γ 2 φ2 1+β21n * 1j s 2a    .
(S28)

The B 3 matrix corresponds to the projection matrix of species 2, which implies that the exclusion equilibrium is not stable when species 2 has a growth rate larger than 1 (i.e., the modulus of B 3 's leading eigenvalue is larger than 1). The stability of this two-dimensional system can be investigated thanks to the Jury conditions,

given by the equations (S29, S30, S31):

1 -tr J(n) + det J(n) > 0 (S29) 1 + tr J(n) + det J(n) > 0 (S30) 1 -det J(n) > 0 (S31)
with J(n) the Jacobian matrix of the two-compartment system whose stability is being evaluated. If these three conditions hold, the system is stable. If the first condition (equation (S29)) is violated, one of the eigenvalues of J(n) is larger than 1, which means in our context that the excluded species can invade the community. The violations of the conditions (S30) and (S31) correspond to the creation of limit cycles [START_REF] Neubert | Density-dependent vital rates and their population dynamic consequences[END_REF]. We then check these conditions on the B 3 matrix. We have:

tr

B 3 = (1 -γ 2 ) φ 2 1 + β 21 n * 1j + s 2a det B 3 = (1 -γ 2 ) φ 2 1 + β 21 n * 1j s 2a -γ 2 φ 2 1 + β 21 n * 1j π 2 1 + α 21 n * 1a .
We start with the third condition (S31):

1 -(1 -γ 2 ) φ 2 1 + β 21 n * 1j s 2a + γ 2 φ 2 1 + β 21 n * 1j π 2 1 + α 21 n * 1a > 0. (S32)
This equation will always be satisfied with biologically meaningful parameters and strictly competitive interaction between species:

• The third term of (S32) is positive because each parameter is positive,

• (1 -γ 2 ) φ2 1+β21n * 1j s 2a < 1 because 0 ≤ γ ≤ 1, 0 ≤ φ ≤ 1 and all β's are positive.

Since tr B 3 > 0, the first Jury's condition implies the second ((S29) ⇒ (S30)). Then, we only have to check the first Jury condition to describe the stability of the system:

1 -(1 -γ 2 )

φ 2 1 + β 21 n * 1j + s 2a + (1 -γ 2 ) φ 2 1 + β 21 n * 1j s 2a -γ 2 φ 2 1 + β 21 n * 1j π 2 1 + α 21 n * 1a > 0 (1 -γ 2 ) φ 2 1 + β 21 n * 1j -(1 -γ 2 ) φ 2 1 + β 21 n * 1j s 2a + γ 2 φ 2 1 + β 21 n * 1j π 2 1 + α 21 n * 1a < 1 -s 2a (1 -γ 2 )φ 2 (1 -s 2a ) < (1 -s 2a )(1 + β 21 n * 1j ) - γ 2 φ 2 π 2 1 + α 21 n * 1a (1 -γ 2 )φ 2 1 1 + β 21 n * 1j + γ 2 φ 2 π 2 1 -s 2a 1 (1 + β 21 n * 1j )(1 + α 21 n * 1a ) < 1.
We finally have the following condition for the stability of the equilibrium:

C 2 1 + β 21 n * 1j + D 2 (1 + β 21 n * 1j )(1 + α 21 n * 1a ) < 1 (S33)
with C 2 = (1 -γ 2 )φ 2 and D 2 = π2γ2φ2 1-s2a .

We obtain the coexistence outcome domains given in Fig. S1 with the β coefficients given in Table S1.

The results are qualitatively similar to those presented in the main text with adults → adults + juveniles. Why is this sub-case of the "metamorphosis life cycle" model more amenable to analytical insights? Supplement D Simulations of metamorphosis model structure (model

2) extended to a S-species context

The species richness at the end of simulations for the extended model 2 with S species are given in Fig. S5.

species, low variance

Number of species at equilibrium We find here that with a slightly different competition structure (adults only compete with adults and juveniles with juveniles), results are only slightly different from those of model 1. They go in the same general direction: it is difficult to obtain a positive effect of opposite competitive hierarchies on species richness at equilibrium. One could have believed that there was something special with [START_REF] Moll | Competition and coexistence with multiple life-history stages[END_REF]'s life cycle structure with separated adults from juveniles, which would have made many-species coexistence despite competition a little bit easier. We found the opposite: for this particular set of simulations, we have not been able to observe more coexistence with opposite competitive hierarchies on α and β, even for very similar

Figure 1 :

 1 Figure1: Outcomes of the model predicted by invasion analysis depending on the coefficients of the inter-specific competition on fertility (α 12 and α 21 ) with β coefficients indicating coexistence (C), exclusion of species 1 (D) and a priority effect (E). The dotted lines correspond to the outcome boundaries of panel A which is used as a reference. The three asterisks stand for the locations of the three parameter sets corresponding to contrasted coexistence scenarios.

Figure 2 :

 2 Figure2: Elasticity of the equilibrium population densities to changes in parameters for the 3 parameter sets of Table2.

Figure 3 :

 3 Figure 3: Histograms of the number of species persisting at equilibrium for 200 permutations of inter-specific competition coefficients. The bar in gray indicates the number of species persisting at time t = 3000 time steps for the original α ij and β ij with Corr(α ij , β ij ) < 0. From left to right, columns correspond the communities with initially 5, 10 and 40 species. From top to bottom, the lines correspond to situations with low (A-B-C), average (D-E-F) and high (G-H-I) variance on the vital rates.

a

  11 a 12 a 21 a 22 b 11 b 12 b 21 b

a

  11 a 12 a 21 a 22 b 11 b 12 b 21 b

a

  Figure A1: Parameter values and coexistence regions of the three parameter sets of Table 2 (A -set 1, Bset 2, C -set 3) chosen to promote different scenarios of coexistence.

Figure A2 :

 A2 Figure A2: Population dynamics of a community with initially 40 species and opposite competitive hierarchies on fertility and juvenile survival. All species have very similar vital rates due to a low variance in the distribution used to draw the parameter. Each line corresponds to one species.

Figure A3 :

 A3 Figure A3: Rank of β ij versus rank of α ij under a negative correlation of β ij and α ij of -0.9.

  Figure A4: Frequencies of outcomes predicted by pairwise invasion criteria for a parameter set with 40 species with very close parameters and a negative correlation between α ij and β ij with φ = 0.9σ

  invasion criteria of pairs of species

Figure A5 :

 A5 FigureA5: Frequencies of outcomes predicted by pairwise invasion criteria -for singe-vital-rate-competition models -for a parameter set with 40 species with very close parameters and bimodal distribution of interspecific coefficients drawn to promote situations of opposite exclusion within each pair of species.

  FigureA6: Frequencies of outcomes predicted by pairwise invasion criteria for a parameter set with 40 species, with very close parameters and permuted competition coefficient that promoted before permutation situations of reciprocal exclusion within each pair of species.

Figure A7 :

 A7 FigureA7: Histograms of the number of species persisting at equilibrium, with initially 40 species, and for 200 permutations of inter-specific competition coefficients. The bar in gray indicates the number of species persisting at equilibrium for the original α ij and β ij with Corr(α ij , β ij ) < 0 and reciprocal exclusion for each pair of species.

Figure S1 :

 S1 FigureS1: Outcomes of the model with a metamorphosis structure predicted by invasion analysis, depending on the coefficients of inter-specific competition on fertility (α 12 and α 21 ) with β coefficients indicating coexistence, exclusion of species 1 and priority effect. The dotted lines correspond to the outcome boundaries of panel A which is used as a reference.

  Figure β 12 β 21 C 0.06 0.06 D 0.11 0.05 E 0.12 0.12 Table S1: β competition coefficients used to obtain contrasting situations of coexistence.

Figure S5 :

 S5 FigureS5: Histograms of the number of species persisting at equilibrium for 200 permutations of inter-specific competition coefficients. The bar in gray indicates the number of species persisting at time t = 3000 time steps for the original α ij and β ij with Corr(α ij , β ij ) < 0. From left to right, columns correspond the communities with initially 5, 10 and 40 species. From top to bottom, the lines correspond to situations with high, average and low variance on the vital rates.

Table 1 :

 1 Parameter sets chosen to obtain equal equilibrium densities but corresponding to contrasted scenarios of coexistence. We use intraspecific coefficients α ii

	1a s 2a	α 12	α 21	β 12	β 21

Table 2 :

 2 Invasion criteria for species 1/species 2, evaluated for models where a single vital rate is affected by competition. R α and R β are invasion criteria respectively for models with competition only on fertility and juvenile survival. Invasion criteria whose values are larger than 1 are in bold.

	)

  α 11 β 11 x 2 + α 12 β 12 y 2 + (α 11 β 12 + α 12 β 11 )xy + (β 11 + α 11 -C 1 α 11 )x + (β 12 + α 12 -C 1 α 12 )y + (1 -C 1 -D 1 ) = 0 α 22 β 22 y 2 + α 21 β 21 x 2 + (α 22 β 21 + α 21 β 22 )xy + (β 22 + α 22 -C 2 α 22 )y + (β 21 + α 21 -C 2 α 21

  + β 11 n 1j (t) + β 12 n 2j (t) the numerator and denominator, leads ton 1j (t + 2) = π 1 φ 1 n 1j (t) 1 + β 11 n 1j (t) + β 12 n 2j (t) + α 11 φ 1 (t) = π 1 φ 1 n 1j (t) 1 + (β 11 + α 11 φ 1 )n 1j (t) + β 12 n 2j (t).This equation gives n 1j (t + 2) as a function of n 1j (t) and n 2j (t) and is equivalent to the classical 2-species Leslie-Gower model (except for the doubled timestep), and therefore we know its behaviour (fixed point), asshown by[START_REF] Cushing | Coexistence of competing juvenile-adult structured populations[END_REF]. Because n 1a (t + 1) only depends on n 1j (t) and n 2j (t) we can deduce that it won't change the dynamical behaviour.

	Going to t + 2
		n 1j (t + 2) =	π 1 n 1a (t + 1) 1 + α 11 n 1a (t + 1) + 0 × n 2a (t + 1)	(S37)
	with n 1a (t + 1) =	φ1n1j (t) 1+β11n1j (t)+β12n2j (t) . This leads to
		n 1j (t + 2) =	φ1n1j (t) 1+β11n1j (t)+β12n2j (t) 1+β11n1j (t)+β12n2j (t) 1 + α 11 π 1 φ1n1j (t)	(S38)
	which, multiplying by 1 (S39)

Acknowledgements

We thank Sam Boireaud for an exploratory numerical study of some two-species models, as well as Olivier Gimenez for discussions and Coralie Picoche for comments on the manuscript. We are grateful to three reviewers for their constructive input.

Statements and Declarations

Code accessibility Computer codes written for analyses are available at GitHub: https://github.com/g-bardon/StructuredModels_Coexistence

Contributions GB wrote the computer code, produced the figures, and performed numerical analyses. Analytic formulas were derived by FB and GB. The first draft was written by GB with help from FB, both

Competing interests and funding

The authors have no relevant financial or non-financial interests to disclose. Funding for Gaël Bardon's internship was provided by the French National Research Agency through ANR Democom (ANR-16-CE02-0007) to Olivier Gimenez.

with C 2 = (1 -γ 2 )φ 2 and D 2 = π2γ2φ2 1-s2a .

This last expression (eq. 15) gives us an invasion criteria that is larger than 1 if, conversely, species 2 is able to invade the community when rare.

Note that if we remove the competition on one vital rate by setting all the competition coefficients associated to this vital rate to 0, we re-obtain the invasion criteria given by [START_REF] Fujiwara | Coexistence of competing stage-structured populations[END_REF] for their model where a single vital rate was affected by competition. This ensures some internal coherence to the analytical results.

Appendix 2 Coexistence regions obtained through sensitivity analysis for the three scenarios of coexistence

To have another, complementary representation of the sensitivities, we estimated coexistence regions for the three parameter sets using the definition given by Barabás et al. (2014b). Note that we arbitrarily use the term region in place of domain to separate the regions found with this sensitivity analysis from the domains of coexistence found through invasion analysis. The method consists in using the sensitivities at the equilibrium for each parameter in order to determine the region of parameter space where species both persists despite changes in parameter values. It is calculated by finding the smallest (positive and negative) perturbation that would lead to the extinction of one of the species from the value of sensitivities and densities at equilibrium, and makes the (crude) approximation that sensitivities are not only valid close to the fixed point but also far away from it. The coexistence regions for the three parameter sets are presented in Fig. A1.

Appendix 3 Additional simulations with S species, S > 2

The meta-parameter sets, used in the parameter distributions from which the vital rate parameters for each species have been drawn, are provided in Table A1. The last two columns correspond respectively to the mean of inter-specific coefficients (for both α ij and β ij ) and to the variance of the inter-specific coefficients (with the covariance between α ij and β ij equal to -0.9 × σ 2 ).

For initially 40 species, a small variance of parameters across species, and opposite competitive hierarchies, we obtained the dynamics illustrated in Fig. A2.

Supplement C Dynamics of single-species and two-species models

Here we explore the dynamics of single-and two-species models in four contrasted life-history configurations (identical to those previously considered by [START_REF] Neubert | Density-dependent vital rates and their population dynamic consequences[END_REF] in single-species models where one vital rate has Ricker density-dependence):

• s a > 0, γ < 1 which corresponds to the parameters used for illustration in the main text

• s a > 0, γ = 1: iteroparity with fast maturation

• s a = 0, γ < 1: semelparity and delayed maturation

• s a = 0, γ = 1: semelparity and fast maturation.

For single-species models (which we use when performing invasion analysis at the exclusion equilibria), only the last case (i.e. semelparity and fast maturation) leads to 2-cycles with the remainder of parameters as in the main text (see figures below). This stays true when multiplying the fecundity by 10-a choice motivated by the fact that high fecundities generally enhance cycling. For the metamorphosis model, the semelparity and fast maturation case leads to a fixed point when initial conditions are strictly positive for juveniles and adults. Indeed, [START_REF] Cushing | Coexistence of competing juvenile-adult structured populations[END_REF] covers the semelparous, fast maturation version of metamorphosis model and showed that while the nontrivial equilibrium for the single-species model is stable on the positive cone, it is not so when starting on the boundary with either n j = 0 or n a = 0, which can generate so-called synchronous two-cycles (Theorem 2.1, [START_REF] Cushing | Coexistence of competing juvenile-adult structured populations[END_REF]).

For two-species models, the case s a = 0, γ = 1 generates 2-cycles but s a = 0, γ < 1 can also generate 2-cycles for some parameter sets, and some different kinds of cycles (see below).

Semeparous -fast maturation

Figure S4: Population dynamics for the two-species models 1, where competition is generated by adults, with the three parameter sets built to generate three different scenarios of coexistence. We consider 4 different life history configurations, always with the initial conditions n i (0) = (0, 50) to allow for all accessible 2-cycles.

For the reader interested in dynamical aspects, we recommend to consult [START_REF] Cushing | Coexistence of competing juvenile-adult structured populations[END_REF]. To facilitate the reading, we make a correspondence below between their parameter notations and ours. Their model writes

with a state-space vector (J t , A t , j t , a t ) where one of the species is written in capitals.

species. This only strengthens the generality of our finding that opposite competitive hierarchies are unlikely to be a general explanation to many-species coexistence.


Aurélie Antunes

Denis Servent

Rémy Béroud

Michel De Waard

Evelyne Benoit

Novel antinociceptive peptides interacting with sodium channels: from their identification to their multi-scale pharmacological characterization

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Pain care is a global public health priority and current pain relievers, although effective in many cases, are too often associated with side effects. A promising approach to manage pain and avoid central side effects is to target the dorsal root ganglia (DRG) neurons, located in the periphery. These neurons are wellknown to express various proteins involved in pain transmission, including the 1.7 subtype of voltage-gated sodium channels (Na V 1.7), a major and genetically validated target for chronic pain treatment. Thanks to their high affinity and selectivity, venom peptides of animal origin have been identified as an original source of antinociceptive molecules.

The objective of my PhD is to identify and characterize, by multi-scale studies, a large number of novel antinociceptive venom peptides possessing high selectivity and potency for Na V 1.7. For this purpose, an electrophysiological screening on cells overexpressing the human Na V 1.7 allowed to identify novel hits from a unique and large venom collection of Smartox Biotechnology. The different studies, currently undertaken on some selected peptides, aim (1) to chemically synthesize these peptides (2) to characterize their selectivity profile and safety issue, notably on cell lines overexpressing Na V subtypes (from 1.1 to 1.8) and on native systems such as the mouse cardiomyocytes and neuromuscular system, (3) to specify their efficiency and mechanism of action on mouse DRG neurons, (4) to evaluate in vivo their analgesic properties in vigilant mice (5) to determine their pharmacophore and optimize their functional properties, and (6) to provide evidence for the clinical potential of the best candidates in correcting de novo identified Na V channelopathies.

Abstract

Novel antinociceptive peptides interacting with sodium channels :

from their identification to their multi-scale pharmacological characterization

Aurélie Antunes

Methodology and Results

Based on a collection of numerous and various venoms, this research strategy makes possible to identify and characterize different toxins of therapeutic interest in the context of their development as antinociceptive agents.

The evaluation of functional properties used here with the cyriotoxin-1a as an exemple in DRG neurons using manual patch-clamp and in mice in vivo reveals that the peptide exhibits marked antinociceptive properties, as evidenced by its effects on TTX-sensitive Na V channels (Na V 1,7) and on tactile and heat sensitivity of animals. This allows us to consider a clinical development of the studied toxin,

  1,2 , Denis Servent 1 , Rémy Beroud 2 , Michel De Waard 3 , Evelyne Benoit 1 1 Service d'Ingénierie Moléculaire pour la Santé, ERL CNRS/CEA 9004, CEA de Saclay, Université Paris-Saclay, Gif sur Yvette, France 2 Smartox Biotechnology, Saint-Egrève, France; 3 Institut du Thorax, Inserm UMR 1087 / CNRS UMR 6291, Nantes, France