Eric Imbert

Vincent Devictor

J'ai

Mickael Hedde

Je Remercie

Alain Laura

Fabrice Pierre

Stan Farah

Alain Valentin

S Majumder

mais aussi Vincent Miele et Stéphane Dray du LBBE (Lyon) ainsi que l'équipe du cluster MBB ici à

R E M E R C I E M E N T S

M'accepteriez-vous comme stagiaire pendant cette période ?

Il y a sept ans (!), je contactais Sonia et Vincent pour faire mon stage de Licence, avec la ferme idée d'essayer l'écologie, dont mon expérience à l'époque se limitait à quelques heures de cours se comptant sur les doigts des deux mains. Loin de moi l'idée que cette expérience allait être le démarrage d'une montagne russe pleine d'émerveillements, m'emmenant dans des contrées chaudes ou froides, me faisant travailler sur des bestioles allant de 5 grammes à deux tonnes, ou encore utiliser des outils aussi simples qu'un papier et un crayon, ou aussi complexes qu'un supercalculateur. Pour avoir su susciter l'étincelle initiale, me communiquer leur enthousiasme et m'avoir fait démarrer dans l'écologie scientifique, je leur suis infiniment reconnaissant.

Je remercie en particulier Sonia pour son encadrement hors normes lors de ces trois années de thèse, tant par son expertise scientifique que ses qualités humaines, et pour la formidable liberté qu'elle a pu m'accorder dans mes recherches. Grâce à cela, j'ai eu la chance de bénéficier de conditions de travail excellentes et d'un accompagnement plus que précieux.

Je tiens à remercier également toutes les personnes avec qui j'ai eu le plaisir de travailler à l'ISEM, et notamment tous les membres de l'équipe BioDICée. Au delà de superbes moments passés au labo ou ailleurs, beaucoup d'entre eux ont une contribution très concrète dans ce travail de thèse, que ce soit pour la saisie de donnée (merci Claire !), ou le travail de terrain (merci Angeles, Virginia, Simon, Alain, Claire, Laura, Pierre, Stan, Thomas et Sonia !). Je n'oublie évidemment pas tous les autres et en particulier Vincent Dufour pour son ethousiasme pour mes projets "technologiques". Plus généralement, je tiens à remercier aussi le personnel de l'ISEM pour son efficacité et sa bonne humeur, à commencer par Agnès Mignot (directrice), mais aussi Sylvio, Fadela, Florence et ceux et celles que j'oublie ici pour leur accompagnement administratif.

Au cours de cette thèse, j'ai bénéficié de nombreuses discussions très fructueuses avec de nombreux scientifiques, à commencer par les membres de mon comité de thèse -Thierry Dutoit, Philippe Jarne,

2.2

The spatial EWS and their expected trends . . . . . . . Following global changes, species populations will disappear, grow or move in space. The species composition and the functioning of ecological communities will change. Anticipating these changes has thus naturally taken a central place in ecology [START_REF] Mouquet | Predictive ecology in a changing world[END_REF].

2.3

The task is daunting, however. While species abiotic requirements can be reasonably used to predict where and when a species can disappear [START_REF] Devictor | Birds are tracking climate warming, but not fast enough[END_REF], community scale keeps being the realm of contingency and sudden reorganizations [START_REF] Lawton | Are There General Laws in Ecology?[END_REF]. Lakes can shift from clear to turbid as a nutrient threshold is crossed (Scheffer et al., 2003). The introduction of a predator or an invasive species can suddenly dramatically reorganize a community [START_REF] Ellis | Long-term effects of a trophic cascade in a large lake ecosystem[END_REF][START_REF] Estes | Trophic Downgrading of Planet Earth[END_REF]. Arid ecosystems can turn into deserts as grazing pressure crosses a threshold [START_REF] Kéfi | Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems[END_REF]. Fish stocks can suddenly plummet [START_REF] Vasilakopoulos | Resilience and tipping points of an exploited fish population over six decades[END_REF]. Plant-pollinator interactions can collapse following external disturbances [START_REF] Lever | The sudden collapse of pollinator communities[END_REF][START_REF] Dakos | Critical slowing down as early warning for the onset of collapse in mutualistic communities[END_REF]. These dramatic responses arise because ecological systems are not simple but complex, i.e. composed of many interacting elements. Species in a community do not live in a vacuum, but are bound together by interactions. Because many ecological communities are diverse, the web of interactions between species may be unfathomably intricate and difficult to analyze -echoing Darwin's vision of the ecological community as an "entangled bank". This is very real issue as, to date, the consequences of the loss or arrival of a species in a given community remain extremely hard to predict: many species introductions for pest control had and still have unpredicted negative consequences [START_REF] Carvalheiro | Apparent competition can compromise the safety of highly specific biocontrol agents[END_REF].

Despite their complexity, ecological interactions can exbhibit specific structural patterns, that determine the response of the whole community to disturbance. For example, ecological systems can exhibit positive feedback loops -in which a characteristic is amplified because it feeds back on what causes it [START_REF] Wilson | Positive-feedback Switches in Plant Communities[END_REF].

Where strong feedback loops are present, the effect of an external pertubation, even if small, can be amplified. For example, the positive interactions between plants in arid ecosystems are known to play a great role in their response to perturbations [START_REF] Kéfi | Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems[END_REF]. Plants facilitate each other by growing in tight clusters that raise the availability of water within them [START_REF] Rietkerk | Self-organized patchiness and catastrophic shifts in ecosystems[END_REF], up to a certain point in stress. Past this point, facilitation cannot counterbalance stressors and an increase in stress, even small, can trigger an irreversible desertification process.

Not all ecological communities are wired in such a way as to show dramatic shifts. Other studies have also shown that ecological interactions can have specific structures that confer robustness to species extinctions [START_REF] Lever | The sudden collapse of pollinator communities[END_REF][START_REF] Verdú | The Nested Assembly of Plant Facilitation Networks Prevents Species Extinctions[END_REF][START_REF] Memmott | Tolerance of pollination networks to species extinctions[END_REF] and allow them to persist in time despite external disturbance [START_REF] Stouffer | Compartmentalization increases food-web persistence[END_REF][START_REF] Thebault | Stability of Ecological Communities and the Architecture of Mutualistic and Trophic Networks[END_REF].

Much research has naturally been devoted to identifying where on this spectrum a given ecological community falls i. e. identifying its potential response to perturbations, and measure its resilience -its capacity to recover after a perturbation.

A first approach has focused on the particular case of ecological systems exhibiting "threshold" responses, where a change in environmental conditions, even small, can produce a wide shift in one or more coarse characteristic of the ecological system (e. g. total biomass).

These often arise when positive feedback loops are present in species interactions [START_REF] Kéfi | When can positive interactions cause alternative stable states in ecosystems?[END_REF]. This body of research, to which this thesis contributes, has produced generic methods and indicators of degraded resilience for this specific case.

This approach is appealing by its simplicity, yet only applies to specific cases where a specific pattern in interactions is already known (feedback loops). For many communities, interactions between species are unknown and must be mapped. Since the early 2000's and the rise of network science [START_REF] Barabási | Network Science[END_REF], many tools and frameworks have been created to do so and analyze ever-more complex interaction networks [START_REF] Morales-Castilla | Inferring biotic interactions from proxies[END_REF][START_REF] Van Borkulo | A new method for constructing networks from binary data[END_REF][START_REF] Poisot | The dissimilarity of species interaction networks[END_REF][START_REF] Delmas | Analysing ecological networks of species interactions[END_REF]. Ecological interactions, in the context of a network, are represented as links between nodes (e.g. links representing trophic interactions between species or nodes in a food web). Identifying specific patterns in these networks that are known to affect the resilience of the community could be a way to anticipate the effect of disturbance on a given ecological community.

However, this view has been criticized [START_REF] De Ruiter | Food Web Ecology: Playing Jenga and Beyond[END_REF] as too simplistic, as the link between disturbance and interactions is, in fact, a two-way street: ecological interactions may affect the response of a given community to perturbations, but interactions themselves change with perturbations (Figure 0.1). Indeed, many striking reorganizations of ecological communities entail a rewiring of ecological interactions on top of a change in species abundances. For example, work on plant-pollinator interactions have shown that not only they can determine the community response to disturbance [START_REF] Memmott | Tolerance of pollination networks to species extinctions[END_REF] but also be affected by disturbance themselves (Memmott et al., 2007b). In plant ecology, decades of work have shown the change of interaction strengths between plants along environmental gradients [START_REF] Callaway | Positive Interactions and Interdependence in Plant Communities[END_REF]. Describing the response of communities to disturbance thus requires monitoring species abundances as well as interactions between species [START_REF] Gray | Ecological networks: The missing links in biomonitoring science[END_REF].

A tight, and two-way relationship thus exists between the structure of ecological interactions, and the resilience and response to perturbations of the community (Figure 0.1). The way these are studied is quickly changing as new methods and increasingly-detailed datasets become available to ecologists. This thesis aims at contributing to the quickly-changing approaches to these topics. ,c). However, interactions can also occur between species and their environment (d) as species often modify local abiotic conditions. For example, in the diagram, a positive feedback loop is present between producers and soil quality (litter, nutrients). Taken altogether, these interactions form interaction networks, where nodes are species.

Ecological interactions between species, or between species and their abiotic environment (Figure 0.2), can determine the response 5 of a given ecological community to disturbance. To describe correctly this response, it is required to introduce a set of useful concepts such as stability and resilience.

Stability and resilience in ecological systems

The concept of stability is notorious for having multiple definitions in ecology, which impedes comparison between studies (Table 0.1).

For example, a survey of the literature by [START_REF] Donohue | Navigating the complexity of ecological stability[END_REF] revealed that in model systems, stability is most often considered to be asymptotic stability, i. e. the rate of return to equilibrium of a system subjected to an infinitesimal perturbation of its state. This definition of stability is purely deterministic and does not involve any notion of temporal variability, which is what is most often considered as stability in empirical studies [START_REF] Donohue | Navigating the complexity of ecological stability[END_REF]. However, two ecological systems could have the same asymptotic stability, but different temporal variabilities, for example because they are subject to different pertubation regimes.

This discrepancy in the use of the stability concept is one example among many, as pointed out by the multiple reviews about its use of in ecology [START_REF] Pimm | Pimm -1984 -The complexity and stability of ecosystems[END_REF][START_REF] Ives | Stability and Diversity of Ecosystems[END_REF][START_REF] Borrelli | Selection on stability across ecological scales[END_REF][START_REF] Donohue | Navigating the complexity of ecological stability[END_REF][START_REF] Grimm | Babel, or the ecological stability discussions: An inventory and analysis of terminology and a guide for avoiding confusion[END_REF]. To help reduce misunderstandings across studies, authors have suggested making more precise stability statement (e. g. [START_REF] Grimm | Babel, or the ecological stability discussions: An inventory and analysis of terminology and a guide for avoiding confusion[END_REF], or including multiple aspects of stability at the same time in a given study [START_REF] Donohue | On the dimensionality of ecological stability[END_REF].

Because it is related enough to sometimes be considered a component of stability [START_REF] Donohue | Navigating the complexity of ecological stability[END_REF], the concept of resilience -put simply, the capacity of a given system to recover after a perturbation -has followed a similar path. Authors such as [START_REF] Holling | Engineering resilience versus ecological resilience[END_REF] have pointed out the dual definition of resilience in the literature, separating out engineering resilience -how fast an ecosystem recovers after a perturbation -from ecological resilience -the maximum amplitude of disturbance that a system can absorb while maintaining its structure. More recently, authors have also sought to clarify further the concept of resilience [START_REF] Hodgson | What do you mean, 'resilient'?[END_REF], and in particular define a standardized way to measure it on data [START_REF] Ingrisch | Towards a Comparable Quantification of Resilience[END_REF].

Despite this diversity of definitions, stability and resilience have proven to be useful concepts to describe and understand the response of ecological systems to perturbations. In particular, stability landscapes (or resilience landscapes, [START_REF] Hodgson | What do you mean, 'resilient'?[END_REF] underpin much work aiming at measuring resilience in empirical systems, and describing their response to perturbations. the theoretical developments around the notion of stability. In this framework, a system's state is described by a set of state variables. For an ecological system comprising a set of different species, the natural state variables could be the set of abundances of every species. The value of these variables is subject to changes because of internal properties of the ecological system -e. g. population growth -but also because of external factors -e. g. environmental conditions.

Put more formally, in this framework, an ecological system is often modelled mathematically by a differential equation of this form:

dX dt = f θ (X) (0.1)
where, X is the set of state variables defining the state of a system. The function f θ , which often depends on a set of parameters θ, defines the intrinsic properties of the system. For example, the following seminal overgrazing model [START_REF] May | Thresholds and breakpoints in ecosystems with a multiplicity of stable states[END_REF][START_REF] Noy-Meir | Stability of Grazing Systems: An Application of Predator-Prey Graphs[END_REF] describes vegetation growing under grazing pressure:

dV dt = rV(1 - V K ) -c V 2 V 2 -V 2 0 (0.2)
The vegetation biomass V (a state variable) here defines the state of the system, while r (reproductive rate), K (carrying capacity), V 0 (half saturation of grazing rate) and c (grazing rate) are parameters.

In this framework, ecological systems can display stable states, which are specific values X * of the state variables at which the state of the system does not change over time (in the absence of noise). These stable states depend on the values of the parameters θ used in the function f θ , and for a given set of parameters, multiple stable states can exist. The above overgrazing model displays such multiplicity of stables states along a gradient of changing values of c (Figure 0.3).

This modelled view of an ecological system is purely deterministic and does not include stochasticity, i. e. chance events or fluctuations that can happen in a given system. In practice, ecological systems are never purely deterministic and stochasticity always contributes to the observed state of a system at a given point in time, hence to its apparent stability.

In a model system with a deterministic response f θ , the parameters θ and the form of the function f define how a system will change given in current state. Assuming a simple model where the state is defined by a univariate continuous value X, the response of the system can be visualized by defining a potential U:

dU dX (θ) = - dX dt (0.3)
Tracing the values of U for different values of X yields a "stability landscape", where stable states are represented by the bottom of valleys which are basins of attraction. Intuitively, the temporal evolution of the system can be described by a ball rolling down this landscape, i. e. towards stable states. For a given value of the system parameters (e. g., grazing rate c in Figure 0.4), a system can have multiple states, as manifested by several valleys in its stability landscape.

Using this landscape analogy, the engineering resilience (i. e. the time of return to equilibrium) is captured by the slope of the potential.

The ecological resilience can be interpreted as the height of the basins of attractions, i. e. how much the system needs to be perturbed to switch to a different state. The shape of this landscape, and how it varies as environmental conditions (parameters) change defines the resilience of the system.

In particular, one central challenge in resilience ecology has been to identify where and when multiple stable states (i. e. multiple valleys in the potential) are present in empirical ecosystems. When they are, irreversible shifts can occur as one of the stable states is lost (Figure 0.4), and stochastic perturbations can have long-term consequences by making a system shift from one state to another. For this reason, much research has sought to find methods that link the characteristics of the stability landscape with empirical patterns [START_REF] Scheffer | Catastrophic regime shifts in ecosystems: Linking theory to observation[END_REF].

potential analysis

Following the potential analogy, the state of a system will always tend towards stable states (the ball rolls down into the valleys of the stability landscape). As a result, it is more likely to be observed close to its stable states than away from them. This intuition has been formalized for a system defined by a single state variable by [START_REF] Livina | Potential analysis reveals changing number of climate states during the last 60 kyr[END_REF].

Considering a simple system defined by a single state variable X, repeated observations of the same system will yield a set of values whose empirical density P x can be computed (e. g. through kernel density estimation). If the temporal response dX/dt of a system is subject to white noise of amplitude σ, then the potential U has a oneto-one relationship with the above empirical density through:

U(x) = - σ 2 2 logP x (0.4)
This relationship implies that the stable states of a given system, i. e. the bottom of valleys in the potential U, should correspond to modes in the empirical distribution of observed states. In particular, if a system has two well-separated stable states, then two modes should be observed in the empirical distribution of states [START_REF] Livina | Potential analysis reveals changing number of climate states during the last 60 kyr[END_REF].

This formal link between observed state and potential has had great success in ecology and many studies have used it to interpret multimodality as potential alternative stable states in ecological systems.

For example, [START_REF] Hirota | Global Resilience of Tropical Forest and Savanna to Critical Transitions[END_REF] identified tropical forest and savanna as alternative stable states that can coexist under similar levels of rainfall, based on the presence of multimodality in the distribution of tree cover. Similar studies have been carried out for a variety of ecological systems such as boreal biomes (Scheffer et al., 2012b), fish stocks [START_REF] Vasilakopoulos | Resilience and tipping points of an exploited fish population over six decades[END_REF], arid systems [START_REF] Berdugo | Plant spatial patterns identify alternative ecosystem multifunctionality states in global drylands[END_REF] or small rock pools [START_REF] Buschke | The dynamics of mountain rock pools -Are aquatic and terrestrial habitats alternative stable states?[END_REF].

Yet the observation of multimodality in empirical systems only weakly supports the presence of alternative stable states. The presence of the latter entails the presence of bimodality in data, but the reverse is not necessarily true. Bimodality in data could be underpinned by the presence or absence of an unmeasured factor that has not been taken into account. The presence of multi-modality, however, remains useful as exploratory analysis to identify where potential alternative stable states -i. e. non-linear responses to perturbationsare present. We use this approach in Chapter 1 to identify potentially fragile plant vegetation communities in the subalpine meadows of Sierra Nevada (California, USA).

critical-slowing down and indicators

As a system gets closer to a critical point, the potential landscape flattens. This is particularly apparent in the fold bifurcation (Figure 0.4). In this case, at the critical point, one stable equilibrium disappears and the system may shift to a potentially very different state (i. e. a valley in the potential "fills up" and disappears). In practice, this translates into a phenomenon called critical slowing-down.

As a system slows down, s pecific aspects of the system are expected to arise in its temporal and spatial properties [START_REF] Scheffer | Catastrophic shifts in ecosystems[END_REF]. Namely, an increase in auto-correlation is expected to occur, because its state at a given time (or position) will tend to be more similar to its state in the previous state (or nearby position). Small perturbations will also have a larger effect as the system takes more time to return to its equilibrium (decreased engineering resilience), producing an increase in temporal and spatial variance. These two metrics were the first suggested in the literature, but others have then been developed since (e. g. skewness [START_REF] Guttal | Changing skewness: An early warning signal of regime shifts in ecosystems[END_REF], informationbased metrics Dakos & Soler-Toscano (2017), etc.).

Measuring these indicators (or early-warning signals) of critical slowing-down on empirical data could in principle inform on when and where a given ecosystem has a decreased resilience, in particular when it happens because of an upcoming shift to another stable state. For this reason, a large body of literature investigates empirical data for such changes in critical slowing down. While they have been observed with reasonable certainty in experiments [START_REF] Veraart | Recovery rates reflect distance to a tipping point in a living system[END_REF][START_REF] Drake | Early warning signals of extinction in deteriorating environments[END_REF][START_REF] Clements | Including trait-based early warning signals helps predict population collapse[END_REF], their relevance for empirical data remains controversed. Studies have found empirical evidence showing an increase in these metrics before observed shifts in temporal data [START_REF] Dakos | Slowing down as an early warning signal for abrupt climate change[END_REF] and spatial data (e. g. [START_REF] Eby | Alternative stable states and spatial indicators of critical slowing down along a spatial gradient in a savanna ecosystem[END_REF]. Yet many of these cases identify shifts after the fact, which questions whether indicators based on them can serve as accurate forecasting tools (Boettiger & Hastings, 2013). Progress on these issues could come from assessing the predictive capacity of these indicators on more empirical datasets. Chapter 2 contributes to these questions by implementing the computation of spatial indicators in common statistical software (R). A more extensive review and summary of indicators is also available there.

Because they are only based on the generic properties of criticalslowing down, the above indicators can potentially apply to a large variety of ecological systems. For this reason, they have received much exposure in the ecological literature. However, other indicators of resilience are based on spatial patterns that arise in specific systems, where interactions between species in a given ecosystem have a specific effect on spatial structure.

specific indicators

One of the main sources of the presence of alternative stable states is the presence of positive feedbacks between elements of a given system. For example, in some arid systems, vegetation tends to favor the local recruitement of new individuals (see also Section 0.2.2.2). In other words, vegetation favors its own recruitment, hence the positive feedback. This has been suggested as a major mechanisms underlying the dynamics of arid systems [START_REF] Kéfi | Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems[END_REF], which are thought to show irreversible transitions associated with the presence of alternative stable states. Because this positive feedback between plant individuals at local scale has a spatial component, it produces a spatial distribution of vegetation with specific characteristics. Some of these characteristics have been shown to vary as stress increases and when an irreversible transition is approaching.

For example, the shape of the distribution of vegetation patch sizes has been suggested as a possible indicator of upcoming irreversible transition, hence decreased resilience (see Chapter 2 and [START_REF] Kéfi | Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems[END_REF][START_REF] Kéfi | More than a meal. . . integrating non-feeding interactions into food webs[END_REF]. Others have also suggested to use the shape of vegetation patches as possible indicators of degradation [START_REF] Rietkerk | Self-organized patchiness and catastrophic shifts in ecosystems[END_REF].

These metrics and indicators try to detect upcoming shifts in an ecological system, for the cases where a particular pattern is present in interactions: a feedback loop. This is an example of a structural pattern in an ecological network that determines how the community will respond to perturbations. However, for most ecosystems, the interactions themselves -between species, or between species and their local environment -are not known. It is thus often necessary to map them in empirical systems.

Mapping interactions in ecological systems

Inferring interactions in ecological communities

The study of interactions between organisms has always been at the core of ecology, but it is only in the past several decades that research has sought to map them systematically in empirical communitiesi. e. drawing a list of all pairwise interactions between species in a given community. An advantage of doing so is that these interactions can be viewed as elements of an interaction network, where nodes represent species and links interactions. By viewing them as networks, a great number of visualisation and statistical tools become available.

Trophic interactions have been the first to be mapped as part of a network [START_REF] Proulx | Network thinking in ecology and evolution[END_REF]. This has changed in the past couple of decades as ecologists have also focused on other types of interactions such as interactions between plants and pollinators or between hosts and parasites [START_REF] Ings | Ecological networks -beyond food webs[END_REF]. Recent studies have also used networks including several types of interactions [START_REF] Kéfi | Network structure beyond food webs: Mapping non-trophic and trophic interactions on Chilean rocky shores[END_REF][START_REF] Pocock | The Robustness and Restoration of a Network of Ecological Networks[END_REF].

Mapping interactions can be straightforward when direct observation is possible (although measuring an interaction strength is often more difficult). For example, plant-pollinator interactions can be mapped by direct observation of pollinators' visits to plants. When direct observation is impossible, using proxies is required to reveal where an interaction is present or not [START_REF] Morales-Castilla | Inferring biotic interactions from proxies[END_REF].

Such a widely-used proxy is the spatial distribution of species. The intuitive idea is that two species that interact in a negative way (e. g.

through competition or by producing hybrids with reduced fitness) should be segregated in space. Conversely, two species that interact positively should often co-occur in the same site.

This approach stems from investigations of the effect of competition on the biogeographical distribution of species. Cody andDiamond's work in 1975 (Cody &Diamond, 1975) was not the first to apply co-occurrence analyses 1 , but marked the beginning of a heated debate on their computation and interpretation in empirical data [START_REF] Sanderson | Patterns in Nature: The Analysis of Species Co-Occurrences[END_REF]. The authors interpreted "checkerboard distributions", where two species rarely ever co-occur on a set of islands, as evidence for competition. This interpretation was disputed by [START_REF] Connor | The Assembly of Species Communities: Chance or Competition?[END_REF] later rejoined by Diamond [START_REF] Diamond | Examination of the "null" model of Connor and Simberloff for species co-occurrences on islands[END_REF]. This debate has spawned along the way a great body of research investigating the link between the geographical distribution of species and their interactions, much beyond the original, species-on-islands problem.

In particular, [START_REF] Connor | The Assembly of Species Communities: Chance or Competition?[END_REF] pioneered the use of null models. Null models are "pattern-generating models" [START_REF] Gotelli | Null Models in Ecology[END_REF], that serve as a reference to produce the patterns that would be expected given the absence of process of interest. For example, in the absence of any interspecific interactions, one could consider that individuals of a pair of species would distribute at random in a set of sampled sites. To this day, null models are still widely used to infer potential interactions from the spatial distribution of species (Box 1).

Historical co-occurrence studies mostly focused on singling out species from a regional pool that potentially competed with each other. For example, Cody & Diamond (1975) mostly focused on competition (among guilds of New Guinea pigeons). Recent co-ocurrence approaches are based on the same principles, but differ in at least two key aspects. First, emphasis is put on positive interactions as well as negative interactions. Second, recent studies do not focus on single interactions but rather view them as part of a network with a given structure. This structure can be analyzed and described, and ecologists have borrowed many tools from network science to do so [START_REF] Delmas | Analysing ecological networks of species interactions[END_REF]. For example, both these aspects are present in a recent co-occurrence framework suggested by Morueta-Holme et al.

(2016).

These methods produce sets of positive and negative associations between species, i. e. association networks 2 (Figure 0.5).

sometimes also named more explicitely "spatial association networks", (e. g. [START_REF] Saiz | Evidence of structural balance in spatial ecological networks[END_REF] Box 1: An example null-model-based analysis of species co-occurrences Classical co-occurrence analyses are done on occurrence (presence/absence) data, that is represented as a matrix of ones and zeros:

sp. 1 sp. 2 sp. 3 sp. 4

site 1 0 0 1 1 site 2 1 1 0 1 site 3 0 1 0 1
First, a metric is used to describe how associated two species are. For example, such a metric can simply be the count of cooccurrences when sampling is discrete, a measure called the natural metric by [START_REF] Sanderson | Patterns in Nature: The Analysis of Species Co-Occurrences[END_REF]. Other metrics with slightly different properties are also used in the literature (e. g. the checkerboard score [START_REF] Stone | The checkerboard score and species distributions[END_REF], χ 2 -statistic (Soliveres & Maestre, 2014), (partial) correlation coefficient [START_REF] Morueta-Holme | A network approach for inferring species associations from co-occurrence data[END_REF]). The metric can be most often computed pairwise between each species independently, yielding observed pairwise associations A ij . Here, using the number of co-occurrences as our association metric, we get the following A ij for the data above:

sp.

1 sp. 2 sp. 3 sp.4 sp. 1 . 1 1 1 sp. 2 1 . 0 2 sp. 3 0 0 . 1 sp. 4 1 2 1 .
Note that for occurrence data, intra-specific co-occurrence cannot be computed, as this requires recording co-occurrence between individuals of the same species (see Chapter 3).

The observed pairwise co-occurrence metric is then compared to a null model. The goal of the null model is to produce a set of null site x species matrices that correspond to the expected spatial distribution of species in the absence of interactions. To do so, many studies use constrained permutations of the original site x species matrix to obtain expected null occurrences of species, but recent approaches also use independent data (e. g. from species distribution models, [START_REF] Morueta-Holme | A network approach for inferring species associations from co-occurrence data[END_REF]).

Regardless of the chosen method, a set of N presence/absence matrices are obtained, from which a distribution of N pairwise association values E ij,1 . . . E ij,N can be computed. Significant pairwise number of co-occurrences are then retained in the original association values A ij , depending on whether they fall in one of the tails of the null distribution. Because this method is akin to a carrying out a test (with α=0.05 in the above example) for each pair of species, correcting for multiple tests is often done using the e. g. Benjamini-Hochberg correction. A strength of association is sometimes also defined as the SES (Standardized Effect Size) of the observed value compared to the null distribution. This process produces an association matrix which holds either significant association strengths (which can be positive or negative), or qualitative values representing a positive or a negative association between each pair of species (Figure 0.5).

Viewing interactions as part of a network has produced many new methods for the inference of interactions from spatial patterns. Many of them entail regressing a model of the whole interaction network instead of relying on separate pairwise analyses. For example, some statistical models aim at modelling the effect of the presence of a given species on the probability of presence of another, on top of possible environmental covariates [START_REF] Pollock | Understanding co-occurrence by modelling species simultaneously with a Joint Species Distribution Model (JSDM)[END_REF][START_REF] Ovaskainen | Modeling species co-occurrence by multivariate logistic regression generates new hypotheses on fungal interactions[END_REF]. This method leads to coefficients potentially reflecting directed species interactions (Figure 0.5). Other model-based methods include Markov networks [START_REF] Harris | Inferring species interactions from co-occurrence data with Markov networks[END_REF][START_REF] Van Borkulo | A new method for constructing networks from binary data[END_REF] which produce undirected networks, or bayesian networks [START_REF] Milns | Revealing ecological networks using Bayesian network inference algorithms[END_REF];

Delalandre & Montesinos-Navarro, 2018) which produces directed networks.

Directed network

Undirected network The same classification applies to association networks derived from co-occurrence analyses, which represent associations as signed, qualitative directed or undirected links (b, e). Weighted networks can be also obtained, in which each association has a given strength (c,f). In most studies, the networks obtained through co-ocurrence analyses are signed and undirected (i. e. of type e or f).

- - + + + b) e) Signed network (qualitative) a) d) Unsigned network +0.3 +0.43 +0.17 -0.67 -0.2 c) -0.1 -0.68 f) +0.17
Despite ever-more refined methods, co-occurrence networks often fail to accurately reflect interaction networks [START_REF] Barner | Fundamental contradictions among observational and experimental estimates of non-trophic species interactions[END_REF][START_REF] Freilich | Species co-occurrence networks: Can they reveal trophic and non-trophic interactions in ecological communities[END_REF]. This is due to a variety of possible reasons:

• The effect of interactions on co-occurrence patterns have been shown to vary with the scale of sampling [START_REF] Mcnickle | Checkerboard score-area relationships reveal spatial scales of plant community structure[END_REF][START_REF] Araujo | The geographic scaling of biotic interactions[END_REF]. Different sampling schemes may thus yield different interaction strengths between species.

• Species often co-occur because they share a given habitat, resulting in positive associations that do not reflect interactions.

When the habitat cannot be considered constant, typically when working at biogeographical scale, controlling for this effect may prove very difficult [START_REF] Morueta-Holme | A network approach for inferring species associations from co-occurrence data[END_REF].

• Other assembly processes such as dispersal limitation may also blur the link between co-occurrence and interactions. Two species may for example never co-occur, simply as a result of their inability to colonize shared sites [START_REF] Connor | The checkered history of checkerboard distributions[END_REF].

Many of these general limits apply to analyses carried out at scales much above the size of the organisms. Working at local scales, where patterns are more strongly driven by biotic processes, and where habitat may be considered constant may yield co-occurrence networks that reflect better ecological interactions between species. This is the main reason why co-occurrence is often used as a proxy for interactions in plant communities.

Plant-plant interaction networks and co-occurrence approaches

Because plants are sessile organisms, the potential effect they have on other species and on their environment tends to be most intense around their canopy. As a result, interactions tend to produce a spatial pattern.

A typical example is the case of "nurse" plants in arid ecosystems.

In these systems, nurse plants are often bushes, which facilitate the recruitement and/or the survival of other species by improving local conditions under their canopy. For example, this improvement can occur as a result of shading (which limits evapo-transpiration [START_REF] Valiente-Banuet | Shade as a Cause of the Association Between the Cactus Neobuxbaumia Tetetzo and the Nurse Plant Mimosa Luisana in the Tehuacan Valley, Mexico[END_REF] or making water from deep soil layers available to plants with shallow roots [START_REF] Horton | Hydraulic lift: A potentially important ecosystem process[END_REF]. In this case, facilitation produces a pattern of spatial association between the nurse plant and their facilitated species (their protégées), because the latter have an increased performance and survival under the former's canopy. It seems thus natural to adopt the inverse perspective and use spatial associations between plants to infer facilitative interactions 3 .

It is noteworthy that not all facilitation is the result of such an assymetric relationship between a nurse plant and facilitated species. For example, improvement in

The accuracy of such co-occurrence methods to infer interactions between species have been challenged in the general case [START_REF] Barner | Fundamental contradictions among observational and experimental estimates of non-trophic species interactions[END_REF][START_REF] Freilich | Species co-occurrence networks: Can they reveal trophic and non-trophic interactions in ecological communities[END_REF]. Yet the sessile nature of plants, along with experimental evidence of correlation between co-occurrence and facilitation [START_REF] Graff | Testing the role of biotic stress in the stress gradient hypothesis. Processes and patterns in arid rangelands[END_REF], suggests that they may be informative for the specific case of plant-plant interactions in the field. This random expectation is derived from permutations of the original dataset [START_REF] López | Co-occurrence Patterns along a Regional Aridity Gradient of the Subtropical Andes Do Not Support Stress Gradient Hypotheses[END_REF] or an analytical approximation [START_REF] Saiz | Changes in semi-arid plant species associations along a livestock grazing gradient[END_REF], 2011;[START_REF] Saiz | Evidence of structural balance in spatial ecological networks[END_REF][START_REF] Alados | Plant-plant interactions scale up to produce vegetation spatial patterns: The influence of long-and short-term process[END_REF].

This method yields a network with positive and negative links between species (nodes), which contrasts with [START_REF] Verdú | The Nested Assembly of Plant Facilitation Networks Prevents Species Extinctions[END_REF] where a bipartite network was used. The network structure can then be described, and its changes along space or in time analyzed.

pollination through a shared pollinator can result in facilitation and spatial association between two plants with similar life forms [START_REF] Johnson | Pollination success in a deceptive orchid is enhanced by co-occurring rewarding magnet plants[END_REF]. We can first consider node-level -i. e., species-level -properties.

Species may have different patterns of interactions within a network.

Species can have variable degrees (numbers of links), and the strength of these links may vary from species to species in weighted networks.

In co-occurrence networks, the degree and strength of the interactions of a given species with others may characterize its role within a community. For example, a nurse species within a network would have many positive links -because it is spatially aggregated with many different species -, i. e. a high degree (high number of positive links).

A high degree here is correlated with a greater influence on the community dynamics. The role of species structuring communities has been explored using such network descriptors (e. g., Stipa tenacissima in arid systems, [START_REF] Saiz | Structure and spatial self-organization of semi-arid communities through plant-plant co-occurrence networks[END_REF].

Intermediate-scale properties can also be considered in a network.

For example, small subgraphs, or motifs, describing interactions between triplets of species have been used to describe the local structure of interaction networks. The relative frequency of these different motifs have been shown to determine network-level properties of food-webs [START_REF] Klaise | The origin of motif families in food webs[END_REF]. In signed, undirected networks such as association networks, a similar, triplet-based description has been used to compute the structural balance of association networks.

Signed networks are expected to be balanced, i. e. have an excess of specific triplets, which has been supported by recent evidence in plant-plant association networks [START_REF] Saiz | The structure of plant spatial association networks is linked to plant diversity in global drylands[END_REF]. Feedback loops are also a type of intermediate pattern which may have a strong influence on the dynamics of the community.

Finally, many network-level characteristics have been shown to deviate from random expectation in empirical networks. In particular, many empirical interaction networks have characteristics that could increase the robustness of the ecological community -i. e. its ability to maintain a high species diversity despite disturbance. In particular, mutualistic networks such as facilitation networks between plants exhibit a higher amount of nestedness (how much specialist species interact with a subset of those with whom generalist species interact), which confers the whole community a higher robustness against random species extinctions [START_REF] Verdú | The Nested Assembly of Plant Facilitation Networks Prevents Species Extinctions[END_REF][START_REF] Bascompte | The nested assembly of plant-animal mutualistic networks[END_REF]. In food webs, modularity (how well a network can be separated in different groups) is thought to have a similar stabilizing effect [START_REF] Thebault | Stability of Ecological Communities and the Architecture of Mutualistic and Trophic Networks[END_REF][START_REF] Stouffer | Compartmentalization increases food-web persistence[END_REF]. High nestedness have been found in plant facilitation networks [START_REF] Verdú | The Nested Assembly of Plant Facilitation Networks Prevents Species Extinctions[END_REF], and it is thought to increase the robustness of the community, but (to our knowledge) no other study explicitly link network-level patterns in plant-plant interactions with the stability of the community.

It is noteworthy that new ways of analyzing and describing networks appear at a quick pace, as evidenced by the rapid increase in methodological papers on network analyses in the last decade [START_REF] Pellissier | Comparing species interaction networks along environmental gradients: Networks along environmental gradients[END_REF][START_REF] Delmas | Analysing ecological networks of species interactions[END_REF][START_REF] Wells | Species interactions: Estimating per-individual interaction strength and covariates before simplifying data into per-species ecological 5 networks[END_REF][START_REF] Graham | Towards a predictive model of species interaction beta diversity[END_REF][START_REF] Poisot | The dissimilarity of species interaction networks[END_REF]. New structural patterns that are relevant for ecology may thus be found in the future as methods progress.

Regardless of the scale they describe in a network, many of these metrics are related to the potential response to disturbance of ecological communities (Table 0.2). However, in most cases, the relationship between networks structural patterns and stability is only supported by theoretical studies (an exception is the study by [START_REF] Gilarranz | Effects of network modularity on the spread of perturbation impact in experimental metapopulations[END_REF], which shows experimentally the effect of modularity on robustness).

This view of interactions as fixed entities determining the response to perturbations has received much criticism (e. g., [START_REF] De Ruiter | Food Web Ecology: Playing Jenga and Beyond[END_REF].

Interactions themselves are known to vary in time, e. g. as a result of annual cycles [START_REF] Olesen | Temporal dynamics in a pollination network[END_REF], but also as the result of perturbations. It is thus necessary to understand how they change after disturbance.

Changes of networks along gradients of stress

Current global changes will impact the abundance of extant species as well as their interactions (Memmott et al., 2007b;[START_REF] Tylianakis | Global change and species interactions in terrestrial ecosystems[END_REF]. To understand how ecological networks may respond to these changes, it is necessary to know what drives their structure across space and time. To this end, much work have thus focused on documenting how interactions vary along environmental gradients [START_REF] Tylianakis | Ecological networks across environmental gradients[END_REF], and developing statistical methods to do so [START_REF] Poisot | The dissimilarity of species interaction networks[END_REF][START_REF] Pellissier | Comparing species interaction networks along environmental gradients: Networks along environmental gradients[END_REF]. Plant ecology has a long tradition of studying how interactions vary along gradients of stress, which predates by a wide margin the use of ecological networks. In plant communities, one of the outstanding debate surrounds the prevalence of competition vs. facilitation in plant communities [START_REF] Callaway | Positive interactions among alpine plants increase with stress[END_REF]. Whether facilitation or competition dominates in a given community is thought to strongly determine the temporal dynamics of a community, by producing potential sharp transitions [START_REF] Kéfi | More than a meal. . . integrating non-feeding interactions into food webs[END_REF] or dampening the effect of stressors [START_REF] Callaway | Positive Interactions and Interdependence in Plant Communities[END_REF][START_REF] Wilby | Shrubs, granivores and annual plant community stability in an arid ecosystem[END_REF]. This balance could thus have a central role in determining how the community responds to external disturbance, and many studies have sought to describe and predict how it varies along gradients of stress [START_REF] Callaway | Positive interactions among alpine plants increase with stress[END_REF][START_REF] Maestre | Refining the stressgradient hypothesis for competition and facilitation in plant communities[END_REF][START_REF] Smit | Inclusion of biotic stress (consumer pressure) alters predictions from the stress gradient hypothesis[END_REF][START_REF] Holmgren | Strong facilitation in mild environments: The stress gradient hypothesis revisited: Strong facilitation in mild environments[END_REF].

With the use of co-occurrence networks in plant ecology, this topic has benefited from a new, correlative approach to quantify the facilitation/competition balance in plant communities. Co-occurrence networks in this context can be summarized into a ratio quantifying the relative frequency of positive and negative links (Saiz & Alados, 2012, , Chapter 3). A positive ratio is most often interpreted as a dominance of processes linked to facilitation over competition in the local community. This new approach, because it does not depend on tedious pairwise experiments to infer plant-plant interactions, could allow documenting a much increased variety of plant communities along environmental gradients. It has been successfully carried out to show that facilitation promotes diversity in arid systems [START_REF] Saiz | The structure of plant spatial association networks is linked to plant diversity in global drylands[END_REF], and that the competition/facilitation balance at the network scale varies along gradients of grazing [START_REF] Saiz | Changes in semi-arid plant species associations along a livestock grazing gradient[END_REF] and aridity [START_REF] Alados | Plant-plant interactions scale up to produce vegetation spatial patterns: The influence of long-and short-term process[END_REF].

Many of these studies based on co-occurrence networks focus on arid ecosystems, where facilitation is known to have a strong effect on community dynamics. Few studies concern more mesic systems, despite an ongoing discussion as to whether facilitation should also be important in those systems [START_REF] Holmgren | Strong facilitation in mild environments: The stress gradient hypothesis revisited: Strong facilitation in mild environments[END_REF]. Using an approach based on co-occurrence networks, Chapter 3 contributes to this debate by focusing on the relative importance of competition vs.

facilitation in plant communities along grazing gradients of La Crau (France).

Perturbations and temporal dynamics in networks

Few studies have focused on the temporal variations of interaction networks. Most of this work has focused on plant-pollinator networks, probably because of their ease of sampling when compared to other types of networks, which makes them easily replicable in time. These studies have shown that these networks do show temporal variations in time throughout the year [START_REF] Caradonna | Interaction rewiring and the rapid turnover of plantpollinator networks[END_REF][START_REF] Olesen | Temporal dynamics in a pollination network[END_REF], and respond favorably to restoration efforts [START_REF] Kaiser-Bunbury | Ecosystem restoration strengthens pollination network resilience and function[END_REF]. In vegetation studies, despite the known existence of year-to-year variations in interaction strengths [START_REF] Tielbörger | Temporal environmental variation tips the balance between facilitation and interference in desert plants[END_REF], no study (to our knowledge) focuses on the temporal variations of interaction networks in plant communities. This is the object of study of Chapter 4, which focuses on the recovery of potential interaction networks in plant communities after disturbance.

Despite new methods [START_REF] Blonder | Temporal dynamics and network analysis: Temporal dynamics and network analysis[END_REF], understanding how ecological networks vary in time, and how they recover from perturbations appears as a new challenge for network-based analyses of ecological communities.

objectives and research questions

As we have seen through this introduction, interaction networks in ecological communities present specific structural patterns (feedback loops, nestedness, etc.) which determine their response to disturbance.

However, disturbance itself may affect the abundance of species and the strength of interactions. This thesis focuses on different aspects of this loop, using plant communities as model systems.

The first chapters focus on a structural pattern which is specifically relevant to ecosystem management and conservation: feedback loops.

When these are present in a community, possible large shifts may occur despite small perturbations, leading to ranges of environmental conditions where ecosystem resilience can be impaired. The two first chapters focus on detecting the effect of a particular structural pattern -feedback loops -on the response to disturbance.

The two remaining chapters adopt the inverse perspective, and evaluate the effect of perturbations on the structure of interaction networks.

In particular, they focus on the effect of perturbation on a specific structural pattern that has been the object of a decade-long debate in vegetation studies: the positive/negative interaction balance. These two chapters are based on an extensive dataset that is the result of field work carried out -in collaboration with many assistants! -over the course of this PhD program in La Crau, France.

Chapter 3 specifically aims at answering the following research question: How does biotic stress alter community-scale positive/negative interaction balance ?. We answer this question using a recent approach based on association networks, and apply it to a dataset documenting plant communities in 16 grazed sites in La Crau, France. Using a chronosequence approach, a temporal perspective is adopted in Chapter 4 to ask How resilient are plant-plant association networks ? How does it compare to other facets of the community ?.

note to the reader

This thesis comprises an introduction, four chapters, and multiple appendices. Chapter 1 and 2 can be read independently, but it is probably best to read Chapter 3 before Chapter 4 as the latter builds upon methods explained in detail in the former. Appendices complement individual chapters with background information about the study sites (Appendix 2), or more extensive mathematical details about specific topics (Appendix 3). Before each chapter, a small transition reminds the reader of where the work is situated in the overarching topic of the link between resilience and ecological interactions in ecological communities, along with some details about the ecological system under focus.

Part II C H

A P T E R S 1 I D E N T I F Y I N G P O S I T I V E F E E D B A C K S I N V E G E TAT I O N C O M M U N I T I E S O F R E M O T E A R E A S : I N S I G H T F R O M S U B A L P I N E M E A D O W V E G E TAT I O N C O M M U N I T I E S context of the work
It is essential for conservation to be able to anticipate how a given ecological system will respond to pertubations. Where ecological systems are resilient, small perturbations in species densities and environmental conditions produce small changes in the state of the system. However, some systems can exhibit large changes in their state, even after a small disturbance. This often occurs where a specific structural pattern exists in the interactions of a given system: feedback loops. It can happen, for example, when two species enhance each other's recruitment (e. g. through facilitation, [START_REF] Kéfi | Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems[END_REF].

Identifying feedback loops from ecological patterns may be challenging. Previous approaches have attempted at doing so based on the criterion of multi-modality: where feedback loops are active, the distribution of observations of a given ecological system may exhibit multiple modes [START_REF] Livina | Potential analysis reveals changing number of climate states during the last 60 kyr[END_REF]. This approach has been carried out on simple, aggregated data (e. g. total forest cover, Hirota et al. 2011), but it remains unclear how it could be extended to multi-species assemblages. We suggest a way to do so in this chapter, and apply it to iconic conservation areas: the national parks of Sierra-Nevada, California.

The Sierra-Nevada is located in inner California (USA) and represents a the dominant mountain range of the state. In these areas, covered mostly by forests of mixed pine species, meadows occur locally wherever topology favors the local accumulation of water [START_REF] Weixelmann | A field key to meadow hydrogeomorphic types for the Sierra Nevada and southern Cascade Ranges in California[END_REF]. These meadows, despite representing a small fraction of the total protected areas, represent islands of essential habitat with specific, endemic species of amphibians such as the Yosemite toad (Anaxyrus canorus) and the Yellow-legged Frog (Rana muscosa).

Many accessible meadows are used for recreation, and others, including remote ones are used to establish grazing camps. This use has been sparked a debate about the appropriate use of meadows in wilderness areas, and called for an assessment of the impact of grazing on subalpine meadows vegetation communities [START_REF] Lee | A multi-scale evaluation of pack stock effects on subalpine meadow plant communities in the Sierra Nevada[END_REF]. Despite the effort to document the impact of short-term stressors on meadow ecosystems, it remains unclear how they will respond to long-term changes such as climate change. It is also unclear as to whether stressors may have a large or irreversible effect on the species composition of meadows.

To provide hypotheses on how these vegetation communities could respond to stressors, we reused and extended previous approaches that aim at identifying non-linear responses in ecological systems [START_REF] Livina | Potential analysis reveals changing number of climate states during the last 60 kyr[END_REF].

I participated to this survey in spring 2012 as part of my Master's degree. This chapter has been submitted as the article 'Plant community patterns reveal areas of low resilience in subalpine meadows', A. Génin, S. R. Lee, E. L. Berlow, S. M. Ostoja, and S. Kéfi. The work presented here is a ongoing revision prior to resubmission.

introduction

It is essential for pro-active management to identify how ecological systems respond to environmental changes [START_REF] Suding | Threshold models in restoration and conservation: A developing framework[END_REF].

In many ecosystems, small changes in environmental conditions or small perturbations produce small changes in the observed state of the ecosystem. However, in others, the same events can result in large and/or irreversible changes in the ecosystem state. It is thus very desirable to identify where and when a given ecological system has such a decreased ability to recover by itself from perturbations, i. e., a decreased resilience.

In principle, identifying precisely how a given ecosystem responds to perturbations requires documenting its temporal response as external conditions are slowly changed, or perturbation events are performed [START_REF] Bestelmeyer | Spatial perspectives in stateand-transition models: A missing link to land management?: Spatial state-andtransition models[END_REF][START_REF] Petraitis | Multiple Stable States in Natural Ecosystems[END_REF]. Such data could correspond to time-series from either observations or experiments [START_REF] Scheffer | Catastrophic regime shifts in ecosystems: Linking theory to observation[END_REF]. However, in most natural ecosystems, obtaining datasets or setting up experiments at the appropriate scales can be both costly and difficult. Instead, researchers often adopt a 'spacefor-time' substitution approach, where the state of a given ecological system is recorded in multiple locations along an environmental gradient, but at a single or several very close points in time ('snapshot' data).

Many studies have used this type of dataset to identify apparent non-linear responses of ecological systems to changes in environmental drivers [START_REF] Hirota | Global Resilience of Tropical Forest and Savanna to Critical Transitions[END_REF]Scheffer et al., 2012a;[START_REF] Vasilakopoulos | Resilience and tipping points of an exploited fish population over six decades[END_REF][START_REF] Eby | Alternative stable states and spatial indicators of critical slowing down along a spatial gradient in a savanna ecosystem[END_REF] by interpreting the changes in the distribution of state values (e.g. distribution of measured tree cover) along an environmental gradient. These studies assume that, on average, a given system is more often observed close to the states that persist the most in time, i. e. its most stable states, than away from them. As a result, at a given point along a gradient, the observed distribution of states of a given system should be centered around its stable states [START_REF] Livina | Potential analysis reveals changing number of climate states during the last 60 kyr[END_REF]. For example, a system with one stable state is expected to show a unimodal distribution centered on a single average value -its stable state (Box 1.1).

When an ecological system exhibits a linear response along a gradient, this average state changes proportionally to the change environmental conditions along the gradient (panel a1 in Box 1.1, [START_REF] Ratajczak | Abrupt Change in Ecological Systems: Inference and Diagnosis[END_REF]. In contrast, when the response is non-linear, small changes in conditions can result in a disproportionate changes in the system's average state, sometimes resulting in a sharp shift as a threshold in environmental value is crossed (panel b1 in Box 1.1).

Some ecosystems can also exhibit multiple modes in their observed distribution of states at a given point along the gradient, i. e. multiple alternative states (panel c1 in Box 1.1). These may be underpinned by external factors, such as an unmeasured variable, which controls whether the observed value is closer to one state or the other. However, internal factors such as feedbacks may also maintain a given ecosystem in one of these alternative states -they are then effectively alternative stable states [START_REF] Petraitis | Multiple Stable States in Natural Ecosystems[END_REF]. In this case, even a short disturbance event may have long-lasting consequences, by making an ecosystem switch durably from one state to another. Large shifts can also occur despite a small change in drivers' values, if the ecosystem is pushed outside of the range where those multiple states exist. Using this modality criterion, we investigated where apparent nonlinear transitions occurred in vegetation communities of subalpine meadows of Sierra Nevada (California, USA). Many of these meadows are located in two major national parks, Yosemite and Sequoia/King's canyon, and exemplify ecological systems under strong scrutiny, where it is very valuable to anticipate the effect of stressors. However, due to their remoteness, it is often impractical to monitor meadows for temporal changes: using snapshot data to infer the potential effect of stressors is thus very relevant in these ecosystems.

Meadows within the Sierra Nevada occupy a small portion of the total area, but they stand out as highly productive riparian islands within a rugged, semi-arid landscape (Barbour et al., 2007). In those meadows, the species present in vegetation communities respond very closely to changes in the moisture content of the soil [START_REF] Allen-Diaz | Water Table and Plant Species Relationships in Sierra Nevada Meadows[END_REF][START_REF] Weixelmann | A field key to meadow hydrogeomorphic types for the Sierra Nevada and southern Cascade Ranges in California[END_REF], making it the dominant environmental gradient at the meadow scale [START_REF] Lowry | Groundwater controls on vegetation composition and patterning in mountain meadows[END_REF]. This general response to water content can compound locally with other factors to determine the response of local vegetation communities, such as To illustrate the link between multi-modality in the observed states and the type of response of a system, we can consider the seminal model of overgrazing by [START_REF] Noy-Meir | Stability of Grazing Systems: An Application of Predator-Prey Graphs[END_REF]. This model is able to produce various types of responses along a gradient of grazing (C values). Non-linear transitions arise when strong positive feedbacks are present in the system (b1, c1), which are visible in the deterministic response of the model (red line). The black points indicate "snapshot" values, which would be observed by sampling along the grazing gradient different instances of this system (with added noise) at the same point in time.

Using a moving-window approach, we can characterize the modality of the distribution of observations using component analysis (described in Methods). When a distribution tends to be more bimodal, the frequency of counts for two components increases, which arises in the vicinity of non-linear responses (b2, c2). Identifying changes in modality in empirical data could improve the understanding of where along an environmental gradient a non-linear response of the system may occur. In such cases, a small change in conditions can trigger a large change in the system's state (blue arrows). An acute disturbance event, such as a reduction in plant biomass, may also trigger a durable change to an alternative state when discontinuous responses are also present (green arrow).

disturbance by animals [START_REF] Klikoff | Microenvironmental Influence on Vegetational Pattern near Timberline in the Central Sierra Nevada[END_REF][START_REF] Lee | A multi-scale evaluation of pack stock effects on subalpine meadow plant communities in the Sierra Nevada[END_REF], soil coarseness [START_REF] Ankenbauer | The effects of soil organic matter on soil water retention and plant water use in a meadow of the Sierra Nevada, CA: Soil organic matter affects plant water use[END_REF] or nutrient content [START_REF] Ratliff | Nutrients in Carex exserta sod and gravel in Sequoia National Park, California[END_REF].

These interactions have the potential to create non-linear responses of the community along the water gradient. For example, the presence of Carex filifolia sods has been correlated with increased soil nutrients 5 [START_REF] Ratliff | Nutrients in Carex exserta sod and gravel in Sequoia National Park, California[END_REF], suggesting a possible positive feedback between veg-etation and soil quality. Reducing vegetation cover in some communities could impair regrowth for a much longer time than elsewhere in meadows. It is thus essential to identify where these occur and their spatial extent in meadows in order to carry out informed management, in particular in the current context of recurrent disturbance such as droughts [START_REF] Belmecheri | Multi-century evaluation of Sierra Nevada snowpack[END_REF][START_REF] Adams | Annually resolved late Holocene paleohydrology of the southern Sierra Nevada and Tulare Lake, California[END_REF].

Using the example of meadows as a case study, we show how and where non-linear responses of plant communities can be identified in subalpine meadows. Combined with an extensive survey of local surface characteristics, we identify possible drivers underpinning these non-linear responses to guide future experiments and management decisions.

1.2 methods

Analyzing changes in modality along a gradient

Characterizing the modality of a given set of observations can be done by analyzing the distribution of mass of a kernel-estimated density [START_REF] Müller | Excess Mass Estimates and Tests for Multimodality[END_REF]. This approach essentially identifies whether a density tends to be composed of multiple 'lumps' or components (multimodal density) rather than a single component (unimodal density), and has been identified as more robust to outliers compared to previous approaches [START_REF] Silverman | Using Kernel Density Estimates to Investigate Multimodality[END_REF].

To analyze changes along the main environmental gradient, we used a moving window approach. For each set of observations of the dataset falling within the window, we first computed the empirical density of the observations with a gaussian kernel of a given bandwidth σ (fixed for the whole analysis). We then rescaled the density so that its maximum was one, and considered n cut thresholds regularly spaced between zero and one. For each threshold, we considered a "slice" of the density and counted the number of observed disjunct components (Figure 1.4). For a given density, this process yields a set of n cut counts of components. A density that is better-described by a

given number of well-separated components will show this number for a higher number of chosen thresholds. How these counts vary along an environmental gradient can inform about changes in the modality of the distribution of states. The typical example could be a density switching from being unimodal to more bimodal (exhibiting then an increased number of counts with two components) around a certain value of the environmental gradient.

This method can be applied to univariate densities (Figure 1, top;

Box 1), which is what is classically done when the state of the system is defined by a single value (e. g. [START_REF] Hirota | Global Resilience of Tropical Forest and Savanna to Critical Transitions[END_REF]. However, it can also be applied to multivariate densities, obtained through multidimensional kernel density estimation (Figure 1b), which is very relevant to cases where the state of the system is not well-defined by a single value (e. g. total biomass) but several (e. g. abundances of individual species).

n cut should be set to a large value to identify differences in height of the density (we used 512 here, which was enough to detect changes in modality along the gradient). We refer to this process as "component analysis" later in this work. noise. On univariate data (a,c), slicing the density yields a single component for a majority of chosen thresholds. On bimodal data, the same procedure yields two components for most thresholds chosen (b, d). Note that very low thresholds identify spurious components because of outliers in the density, but this effect disappears quickly as higher thresholds are used

Meadow sampling

Sampling in subalpine meadows was carried out as described in [START_REF] Lee | A multi-scale evaluation of pack stock effects on subalpine meadow plant communities in the Sierra Nevada[END_REF]. Within each meadow, surveys were carried out in 2x2m plots, regularly spaced within belt transects laid out every 40m perpendicularly to the main meadow drainage. Ocular estimates of 28 upper soil characteristics were measured (Appendix 2). We measured volumetric moisture content in the soil top 12 cm, as averaged from 5 measurements of a TDR moisture probe. Ocular cover estimates of each plant species was recorded in 8 25x25cm subplots and averaged to obtain an estimation of the relative abundance of each species in the plot. We then analyzed the response of vegetation communities -in terms of species composition -along the natural moisture gradient in meadows. To attenuate the effect of year-to-year variations in precipitation on the measure of the moisture gradient, we standardized the measured water content within year. Hereafter, water content values are given in standard deviations (s.d.).

Empirical data analyses

The state of meadow vegetation communities is defined by the relative abundances of all species. Because carrying kernel density estimation cannot be done on data with many dimensions (typically, computational costs become too high after 5-10 dimensions [START_REF] Blonder | The n-dimensional hypervolume[END_REF], well below the total number of observed species), we had to summarize the variations of vegetation species compositions to a few axes using ordination. In principle, using Principal Coordinate Analysis (also named metric Multi-Dimensional Scaling) is the method of choice as it preserves as much as possible the dissimilarities between sites. However, it is often empirically found to produce strong arch artifacts on vegetation data, in particular when species turnover is high along a gradient [START_REF] Legendre | Numerical Ecology. Number 24 in Developments in environmental modelling[END_REF]. Non-metric dimensional scaling could be considered as an alternative but, because it preserves only the ranks of the original dissimilarities between sites, the resulting axes summarizing species turnover cannot be used quantitatively. We therefore chose a Correspondence-analysis-based method, Detrended-Correspondence Analysis, which, despite criticisms [START_REF] Jackson | Putting Things in Order: The Ups and Downs of Detrended Correspondence Analysis[END_REF], (i) produces axes that can be used as quantitative variables and (ii) is not subject to arch-like artifacts.

We carried out sensitivity analyses to confirm that our results are robust to the use of other dimensionality reduction methods 2 .

We used the two first axes of the DCA and carried out a component analysis in a moving window along the moisture gradient (we set the window size to 0.5). For the component analysis, we used the bivariate diagonal density estimation kernel automatic bandwidth selection (bandwidth.nrd in package MASS). We carried out the component analysis 100 times for each window subset, using a random subsample of 150 of the observations falling within the window to make sure sample size for density estimation was constant along the water gradient.

Because the component analysis identified two ranges of water content with increased support for two components, dry areas and wet areas, we extracted these subsets from the whole dataset and investigated soil attributes that could underpin these bimodal distributions.

We refer to dry areas for sites where water content was below -1 s.d.

and to wet areas for those where water content was above 1 s.d. We clustered each subset into two groups based on species composition, using hierarchical clustering with a Ward distance on a Bray-Curtis distance matrix. We contrasted the environmental covariates of each group to discuss possible underlying mechanisms for the apparent non-linear response along the water-content gradient. We used twosample Mann-Whitney tests (wilcox.test in R) to contrast the distribution of each plot characteristics across groups, adjusting the p-values for multiple tests using the Benjamini-Hochberg correction. Because sampling was regularly-spaced in meadows, we used the frequency of each group as a proxy for its spatial extent in the meadows.

We carried out the analyses using R (version 3.5.0), along with the R package 'vegan' v2.5 [START_REF] Oksanen | Vegan: Community Ecology Package. R package version 2[END_REF]), 'ks' v. 1.11.1 (Duong, 2018) and MASS [START_REF] Venables | Modern Applied Statistics with S. Statistics and computing[END_REF].

results

In a proof-of-concept model (Box 1.1), simulations show that multiple components arise in cases of non-linear responses of the system, and that they are well-identified by our approach. Empirical data anal- It also showed a significant difference in water-regime attributes (cover of dry and moist soil), and disturbance attributes (stock prints, rodent mounds), although they had a very small estimated cover difference (below 10 -4).

discussion

Patterns suggested two ranges of the moisture gradient where multimodality of species composition highlighted apparent non-linear re- sponses of vegetation to the moisture gradient. These two modes in species composition are correlated to differences in local abiotic conditions and animal disturbances. This pattern could be driven entirely by independent abiotic differences, but at this scale, vegetation itself is also likely to influence the local conditions through positive feedbacks [START_REF] Wilson | Positive-feedback Switches in Plant Communities[END_REF]. Disentangling causal drivers of multimodality in these areas would require work beyond observational analyses, but it is important to acknowledge this possibility in a conservation context. If vegetation-environment feedbacks are active, i. e.

the response of the local vegetation is linked to intrinsic, possibly nonlinear dynamics, and not entirely driven by external factors. A small change in annual moisture, or a strong disturbance event could thus lead to unexpected shifts in species composition and different ecosystem functioning. We discuss this possibility in the next paragraphs.

Possible ecological mechanisms behind bimodality in meadows

In wet areas of meadows, our analysis identified a Carex-vesicaria- In dry meadow areas, the analysis identified a Carex filifolia-dominated community contrasting with a more diverse mesic community.

dominated
The C. filifolia community was observed to have coarser upper soil particles (sand, pebble), and a higher amount of bare ground. Significant differences were also found in terms of disturbances by animals and water regime. The presence and characteristics of dry communities in meadows is known to respond to external parameters such as rodent disturbance, erosion processes and soil nutrients [START_REF] Klikoff | Microenvironmental Influence on Vegetational Pattern near Timberline in the Central Sierra Nevada[END_REF].

Plants growing in these areas are known to also influence in turn some of these parameters. For example, the sods formed by C.filifolia are known to interact with rodents [START_REF] Klikoff | Microenvironmental Influence on Vegetational Pattern near Timberline in the Central Sierra Nevada[END_REF], and can improve local conditions by trapping smaller soil particles which hold more nutrients [START_REF] Ratliff | Nutrients in Carex exserta sod and gravel in Sequoia National Park, California[END_REF] and water [START_REF] Wood | Holocene Stratigraphy and Chronology of Mountain Meadows[END_REF][START_REF] Ankenbauer | The effects of soil organic matter on soil water retention and plant water use in a meadow of the Sierra Nevada, CA: Soil organic matter affects plant water use[END_REF]. The possible presence of this feedback means that acute disturbance of these areas could have much longer-termed effects than expected and possibly require restoration. Further work focusing on identifying the type of dynamics at play in these dry areas would yield a better understanding of their resilience to local degradations and climate change.

In both wet and dry areas, water-related parameters stood out as possible drivers of non-linear responses. With the increase of temperatures associated with climate change, snow/rain ratios and total spring snowpack are decreasing across the Sierra Nevada [START_REF] Knowles | Trends in Snowfall versus Rainfall in the Western United States[END_REF][START_REF] Barnett | Human-Induced Changes in the Hydrology of the Western United States[END_REF]. These long-term shifts have the potential to extensively alter meadow hydrology, to which the identified wet and dry communities may be the most sensitive. Short but acute events such as drought could produce persistent shifts of mesic communities, which tend to be the most diverse, into less-diverse dry communities. The Sierra Nevada has recently experienced one of these events [START_REF] Belmecheri | Multi-century evaluation of Sierra Nevada snowpack[END_REF], which included the four driest years (2012-2015) of the last 2000 years [START_REF] Adams | Annually resolved late Holocene paleohydrology of the southern Sierra Nevada and Tulare Lake, California[END_REF]. While the repercussions of the drought have clearly manifested across Sierra Nevada forest ecosystems [START_REF] Potter | Landsat Image Analysis of Tree Mortality in the Southern Sierra Nevada Region of California during the 2013-2015 Drought[END_REF], the impacts on meadow eco-

systems have yet to be thoroughly evaluated. Monitoring changes and setting up experimental research in areas were non-linear responses occur could inform us greatly about such impact.

Method discussion

In this work, we used a descriptive method to reveal where a complex ecological system exhibits apparent non-linear responses, in the absence of temporal monitoring. This method can be readily applied to existing datasets provided that a measure, even noisy, of the main environmental gradient is carried out. To differentiate whether apparent non-linear responses are underpinned by potential positive feedbacks, it remains essential to also have local, detailed data.

We chose to base our component analysis on the two first axes of the ordination because visual inspection revealed the two clusters (Figure 1.5 top), and because we focused on broad variations in species composition (which are contained in the first few axes of the ordination). A similar analysis could be based on a multivariate density estimation that uses more dimensions (in this case, more axes of the ordination). This could in principle reveal more ranges of the environmental gradient with multimodal distributions. However, observations become increasingly distant from each other (in a multivariate space) as dimensions are added. As a result, the possible increased power brought by added dimensions may come with an increased amount of false positives -i. e. ranges of the gradient with increased multimodality, which does not correspond to ecologically relevant differences in the field. Using only one or two dimensions is a conservative option when investigating multi-modality [START_REF] Blonder | Do Hypervolumes Have Holes?[END_REF]. Further work investigating this trade-off, along with computationally efficient methods [START_REF] Blonder | The n-dimensional hypervolume[END_REF] would certainly provide insights to guide the choice of the number of dimensions.

Future refinements of this work should also include sensitivity analyses regarding the choice of ordination, the window sizes and density bandwidths.

Conclusion

In the current context of rapid change and abundant data, a rising number of methods and indicators aim at linking patterns with the stability properties of a given ecosystem [START_REF] Clements | Indicators of transitions in biological systems[END_REF],

based on e.g. time-series properties (Dakos et al., 2014), spatial patterns (Génin et al., 2018a) or in this case, the distribution of states.

Although they do not replace experiment-based ecological knowledge (Boettiger & Hastings, 2013), they can provide valuable monitoring tools to anticipate the response of ecological systems to global changes as well as local disturbances, allowing managers and researchers to make more informed decisions in uncertain contexts. Early-Warning Signals (hereafter EWS) in the literature [START_REF] Scheffer | Early-warning signals for critical transitions[END_REF]. EWS provide a relative measure of the proximity to a critical point, and therefore of a possible ecosystem transition to a different state. Computing and monitoring them could therefore help detecting upcoming critical transitions in ecological systems [START_REF] Scheffer | Early-warning signals for critical transitions[END_REF].

F O S T E R I N G T H E D E T E C T I O N O F E C O S Y S T E M S H I F T S B Y I M P L E M E N T I N G S PAT I A L E A R LY-WA R N I N G S I G N A L S context
Despite their promises, a number of studies have also stressed the weaknesses of EWS. EWS provides information about whether as system is degrading or not, but they do not inform on how far the studied system is from a critical point (but see to be particularly promising because they could allow the systematic monitoring of ecosystems' degradation level based on a few spatial snapshots [START_REF] Kéfi | Early Warning Signals of Ecological Transitions: Methods for Spatial Patterns[END_REF].

To help fuel investigations around spatial EWS, we suggest here a workflow to guide users of a new package for R (R Core Team, 2018), 'spatialwarnings', to facilitate their application to spatial data from ecological and other complex systems. The package can readily compute the spatial indicators from data, assess their significance, and display the results, all in a convenient format much similar to other usual R tasks (e.g. linear regression).

In the following sections, we first give the required background knowledge on spatial indicators. Using case-studies, we then describe how they can be applied to real and model systems with 'spatialwarnings' (version 1.1). For conciseness, we do not show all textual output in the article, but Appendix 2 presents the code used here in more details.

the spatial ews and their expected trends

Current spatial EWS available in the literature come mainly from two different theoretical backgrounds: some are based on the changes of dynamics before a critical point (namely 'critical slowing down'; see below), while others derive from the characteristics of the patchiness of certain ecosystems.

As an ecological system approaches a critical point (Box 2.1), it is expected to take more time to recover from small perturbations [START_REF] Scheffer | Early-warning signals for critical transitions[END_REF], a phenomenon known as Critical Slowing Down (CSD). This is expected to yield an increase in spatial heterogeneity (the spatial variance) of its state variables, as it stays longer away from its average state [START_REF] Guttal | Spatial variance and spatial skewness: Leading indicators of regime shifts in spatial ecological systems[END_REF]. This deviation can sometimes be biased towards higher or lower values of the state variable, resulting in an additional rise of spatial skewness [START_REF] Guttal | Spatial variance and spatial skewness: Leading indicators of regime shifts in spatial ecological systems[END_REF]. For example, a forest close to a transition towards a treeless state could exhibit higher spatial variance in tree cover and a bias towards low cover values (higher skewness) as external stress increases. In addition, because of CSD, a system close to a transition should exhibit higher near-neighbor correlations [START_REF] Dakos | Spatial correlation as leading indicator of catastrophic shifts[END_REF].

These three metrics -spatial variance, skewness and autocorrelation -can, in principle, reflect an upcoming critical point in any dynamical system, hence their reference as 'generic EWS' in the literature [START_REF] Scheffer | Early-warning signals for critical transitions[END_REF].

Later work has produced EWS based on spectral properties (Carpenter & Brock, 2010), hereafter referred to as 'spectrum-based EWS'.

Because of the slowing down, two points at a given distance from each other should be more similar in systems close to critical points [START_REF] Kéfi | Early Warning Signals of Ecological Transitions: Methods for Spatial Patterns[END_REF]. This change in correlation length can be measured by computing the r-spectrum of the data [START_REF] Couteron | Quantifying change in patterned semi-arid vegetation by Fourier analysis of digitized aerial photographs[END_REF]. The r-spectrum is based on the Discrete Fourier Transform and decomposes the spatial data into series of periodic sines and cosines of given wavelengths (the inverse of wave frequency) and amplitude. The amplitude coefficient associated to each wavelength characterizes the dominance of this scale in the spatial data. Close to a critical point, a relative increase in the amplitudes of long wavelengths (lower frequencies) is expected compared to short wavelengths (higher frequencies). To summarize this into a single value, one can use the ratio of the average amplitude of low frequencies over the average amplitude of high frequencies [START_REF] Biggs | Turning back from the brink: Detecting an impending regime shift in time to avert it[END_REF], the Spectral Density Ratio, which is expected to increase along degradation trends.

Additional EWS have been suggested for ecosystems that exhibit a clear spatial structure, such as drylands or savannas [START_REF] Rietkerk | Self-organized patchiness and catastrophic shifts in ecosystems[END_REF]. Characteristics of these spatial patterns, and in particular the distribution of their patch sizes [START_REF] Kéfi | Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems[END_REF][START_REF] Kéfi | Robust scaling in ecosystems and the meltdown of patch size distributions before extinction: Patch size distributions towards extinction[END_REF], have been suggested to be promising indicators of degradation [START_REF] Kéfi | Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems[END_REF][START_REF] Lin | Spatial vegetation patterns as early signs of desertification: A case study of a desert steppe in Inner Mongolia, China[END_REF][START_REF] Berdugo | Plant spatial patterns identify alternative ecosystem multifunctionality states in global drylands[END_REF], hereafter referred to as 'patch-based EWS'. In these ecosystems, when stress is low, density is typically high and there can be a patch whose width or height equals the entire area considered, i.e. a spanning cluster. As density decreases, this spanning cluster breaks down into smaller ones until a barren ecosystem state is reached [START_REF] Kéfi | Robust scaling in ecosystems and the meltdown of patch size distributions before extinction: Patch size distributions towards extinction[END_REF][START_REF] Van Den Berg | Sharpness of the percolation transition in the twodimensional contact process[END_REF][START_REF] Corrado | Early warning signals of desertification transitions in semiarid ecosystems[END_REF][START_REF] Sankaran | Patchiness and scale-free correlations: Characterising criticality in ecosystems[END_REF]. When using spatio-temporal data, this sequential process is reflected in changes of the patch size distribution (PSD) [START_REF] Kéfi | Robust scaling in ecosystems and the meltdown of patch size distributions before extinction: Patch size distributions towards extinction[END_REF]. More precisely, it is expected that, A recent study has further suggested that this gradual truncation of the patch size distribution could be simplified into a single metric [START_REF] Berdugo | Plant spatial patterns identify alternative ecosystem multifunctionality states in global drylands[END_REF]. The idea relies on the fact that many empirical distributions resemble power-laws only for patches above a certain size x min . This power-law range (PLR) can be expressed as

PLR = log(x max ) -log(x min ) log(x max ) -log(x smallest )
where x smallest represents the smallest and x max the largest patch size observed in the dataset. This metric is expected to decrease as the ecosystem becomes degraded.

It is noteworthy that the interpretation of the expected trends in patch-based EWS depend on the underlying ecological mechanisms generating the patterns; these have been particularly well studied in drylands [START_REF] Rietkerk | Self-organized patchiness and catastrophic shifts in ecosystems[END_REF][START_REF] Kéfi | Robust scaling in ecosystems and the meltdown of patch size distributions before extinction: Patch size distributions towards extinction[END_REF][START_REF] Manor | Facilitation, competition, and vegetation patchiness: From scale free distribution to patterns[END_REF].

Therefore, as for all the other EWS, their use without a good understanding of the ecological mechanisms driving spatial structure could lead to misinterpretations of the trends, since altered or even opposite trends can be observed, both in empirical [START_REF] Weerman | Changes in diatom patch-size distribution and degradation in a spatially self-organized intertidal mudflat ecosystem[END_REF] or model systems [START_REF] Pascual | Criticality and disturbance in spatial ecological systems[END_REF][START_REF] Schneider | Spatially heterogeneous pressure raises risk of catastrophic shifts[END_REF].

workflow and ews computation

Because working with spatial EWS can be challenging and errorprone, we provide a workflow for the R package 'spatialwarnings' to make this type of analyses more accessible and reproducible. Our

System state

A behavior of interest in ecological systems are discontinuous transitions via Fold bifurcations [START_REF] May | Thresholds and breakpoints in ecosystems with a multiplicity of stable states[END_REF][START_REF] Noy-Meir | Stability of Grazing Systems: An Application of Predator-Prey Graphs[END_REF], where the prediction of the proximity to a critical point is particularly relevant. In this case, the system of interest has alternative stable states over a range of values of the control parameter (upper and lower solid lines in panel a.1; e.g. high vs low vegetation cover in a rangeland), separated by an unstable equilibrium (dashed line in panel a.1). If a change in the control parameter (e.g. grazing intensity) pushes the system state across the critical point (black dot), the system exhibits a 'catastrophic shift' (black arrow) and collapses from its current state (e.g. high cover) to the other one (e.g. low cover). These transition arise because of the presence of strong positive feedbacks. They can sometimes be present, driving the system close to Foldbifurcation-type dynamics and leading to a non-linear response, but not strong enough to produce a discontinuity (panel a.2) [START_REF] Kéfi | Early warning signals also precede non-catastrophic transitions[END_REF].

Other types of critical points exist without wide shifts in ecosystem states, such as in the case of transcritical bifurcation (panel b), for which the typical example in ecology would be the extinction of a species [START_REF] Clements | Including trait-based early warning signals helps predict population collapse[END_REF][START_REF] Drake | Early warning signals of extinction in deteriorating environments[END_REF].

In all these cases, EWS are expected to reflect the approach of critical points (panels a. The workflow comprises four steps: (i) prepare the data and check for periodicity, (ii) compute the EWS from a list of matrices, (iii) assess In R, 'spatialwarnings' can be installed from CRAN and loaded using the following commands:

install.packages("spatialwarnings") library(spatialwarnings)

Preparation of input data and periodicity checks

Typical raster data, e.g. from airborne or satellite-based remote sensing imagery, often come as image files with multiple bands. These images can be read in R using standard packages (e.g. 'tiff', 'jpeg') to obtain an array (a three-dimensional set of values), or a more complex raster object (package 'raster'). Indicator functions of the package work on matrix objects, the common denominator of all these data types. A conversion of a multi-band raster object to a matrix object is often required before computing EWS, a procedure that depends on the type of data, the system of interest and the questions addressed (see Appendix 3 and e.g. [START_REF] Liu | Image Processing and GIS for Remote Sensing: Techniques and Applications[END_REF] for a more complete reference). In what follows, we assume that the data have been transformed into a matrix object, i.e. a raster surface of a single variable.

As with many spatial analyses, the resolution of the raster data should be such that it appropriately captures the scale of the changes in the spatial structure. For example, if the scale of the spatial structure is small, a coarse raster image may not be able to capture its variations.

In general, using prior knowledge or carrying out preliminary analyses at different scales when the data is available [START_REF] Lam | On the issues of scale, resolution, and fractal analysis in the mapping sciences[END_REF] can help guide the choice of resolution. In addition, the input raster image must be large enough so that the spatial characteristics measured by the indicators are accurate. Subsetting the input data, where possible, could provide a way of estimating measurement errors.

The indicator functions can also operate on lists of matrices to directly obtain trends. This can be the case when one does not only have a single image of the system under study but several ones corresponding to different moments in time or to different instances of a given ecosystems along an environmental gradient. All indicators can be computed on matrices containing binary (boolean, TRUE/FALSE)

data, but some of them can also be computed on continuous (numerical) data (Figure 2.1). A continuous matrix can be transformed into a binary matrix by thresholding or classification. However, this transformation needs to be done carefully depending on the ecological context of the data and is therefore left to the user.

Once a suitable matrix object is obtained, it is necessary to identify whether the spatial structures in the matrix are periodic or not [START_REF] Kéfi | Early Warning Signals of Ecological Transitions: Methods for Spatial Patterns[END_REF]. Some ecosystems, such as patterned bushes in the Sahel, show characteristic changes in spatial periodicity along degradation gradients [START_REF] Rietkerk | Regular pattern formation in real ecosystems[END_REF], which mask the expected increase in EWS before critical points [START_REF] Dakos | Slowing Down in Spatially Patterned Ecosystems at the Brink of Collapse[END_REF]. To identify whether an ecosystem shows periodic patterns, a solution is to detect possible modes in the r-spectrum of the image (Appendix 2, Couteron ( 2002)). The absence of strong modes in the r-spectrum of the input matrices means that the patterns are not periodic and that EWS provided by 'spatialwarnings' can be computed [START_REF] Kéfi | Early Warning Signals of Ecological Transitions: Methods for Spatial Patterns[END_REF]. Specific indicators exist for periodic ecosystems [START_REF] Deblauwe | Environmental modulation of self-organized periodic vegetation patterns in Sudan[END_REF] and may be part of a future 'spatialwarnings' release.

Generic EWS

To illustrate the quantification of generic EWS using the 'spatialwarnings' package, we reproduce here an analysis performed by Eby et al. plot(serengeti.test, along = serengeti.rain) Additionally, individual r-spectra can be displayed for each input matrix using plot_spectrum() (Supplementary figure 1):

• • • • • • • •• • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • •• • • • • Mean
serengeti.test <-indictest(serengeti.spec) plot_spectrum(serengeti.test)
All observed r-spectra should show a near-monotonic decrease: a spectrum with a non-zero mode implies the presence of a periodic pattern in the matrix (Appendix 2), in which case the use of the generic and spectrum-based EWS can be unreliable [START_REF] Dakos | Slowing Down in Spatially Patterned Ecosystems at the Brink of Collapse[END_REF][START_REF] Kéfi | Early Warning Signals of Ecological Transitions: Methods for Spatial Patterns[END_REF], while the shape and size of the patterns may provide more valuable information [START_REF] Kéfi | Early Warning Signals of Ecological Transitions: Methods for Spatial Patterns[END_REF].

Patch-based EWS

We showcase the use of the patch-based EWS on model data from a Forest Gap model by [START_REF] Kubo | Forest spatial dynamics with gap expansion: Total gap area and gap size distribution[END_REF]. In this model of forest spatial dynamics, spatial patterns emerge because of near-neighbor interactions between trees. The dataset, forestgap, contains the model outputs along several values of increasing stress, a list of a binary matrices in which TRUE values identify cells with trees. The data frame forestgap.pars contains the parameters used in the simulations (d, delta and al pha).

EWS can be computed using the patchdistr_sews() high-level function, for example using BIC as criterion to select the best-fitting distribution:

forestgap.patch <-patchdistr_sews(forestgap, best_by = "BIC")

The percolation status, patch distribution type and PLR can be displayed as text using summary(), or as a graph using plot() (Figure 3

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • d: 0.
T H E R E S P O N S E O F H E R B A C E O U S P L A N T -P L A N T A S S O C I AT I O N S A L O N G A B I O T I C S T R E S S G R A D I E N T context
Many interaction networks exhibit patterns that affect the resilience of the whole community. For example, in some plant communities [START_REF] Verdú | The Nested Assembly of Plant Facilitation Networks Prevents Species Extinctions[END_REF], the fact that interaction are nested may prevent species extinctions and maintain a higher diversity. A dominant role of facilitation in a given community has also been identified as the potential source of sharp responses to increase in stress, because they are a necessary condition for positive feedback loops [START_REF] Kéfi | When can positive interactions cause alternative stable states in ecosystems?[END_REF]. Identifying whether the dynamics of a community are more dominated by competition or facilitation could thus inform greatly about the potential response of a given vegetation community to perturbations [START_REF] Callaway | Positive Interactions and Interdependence in Plant Communities[END_REF]. This is one of the reasons

introduction

In many plant communities, much work has shown the fundamental effect of inter-specific interactions on the response of the community to perturbations [START_REF] Verdú | The Nested Assembly of Plant Facilitation Networks Prevents Species Extinctions[END_REF]. This concerns not only competitive interactions, but also positive interactions. For example, in semi-arid communities, nurse effects are prevalent, in which a bush promotes the recruitment of other herbaceaous species under its canopy [START_REF] Valiente-Banuet | Shade as a Cause of the Association Between the Cactus Neobuxbaumia Tetetzo and the Nurse Plant Mimosa Luisana in the Tehuacan Valley, Mexico[END_REF][START_REF] Horton | Hydraulic lift: A potentially important ecosystem process[END_REF].

Such facilitative effects have been suggested as the drivers of shiftlike and irreversible responses of plant communities to perturbations [START_REF] Kéfi | Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems[END_REF]. Knowing where and when facilitative interactions are predominant in a community could inform greatly about its potential response to perturbation.

Species interactions can change along environmental gradients, either directly or as the result of species turnover. For example, a general result seems to be that interactions between plant species tend to be facilitative in abiotically-stressed environments, but switch to competition in less-stressed environments. This result was coined as the Stress-Gradient Hypothesis (SGH) two decades ago [START_REF] Bertness | Positive interactions in communities[END_REF][START_REF] Callaway | Positive interactions among alpine plants increase with stress[END_REF], then further refined by later studies to include different type of stressors and plant strategies [START_REF] Maestre | Refining the stressgradient hypothesis for competition and facilitation in plant communities[END_REF]. For example, the net interaction between a pair of plants can be altered by the presence of herbivores [START_REF] Smit | Inclusion of biotic stress (consumer pressure) alters predictions from the stress gradient hypothesis[END_REF]. Associative protection can arise in grazed conditions, an interaction in which a plant benefits from protection against herbivores due to the repulsive effect of unpalatable neighbors [START_REF] Baraza | Conditional outcomes in plant-herbivore interactions: Neighbours matter[END_REF][START_REF] Graff | Testing the role of biotic stress in the stress gradient hypothesis. Processes and patterns in arid rangelands[END_REF][START_REF] Bagousse-Pinguet | Release from competition and protection determine the outcome of plant interactions along a grazing gradient[END_REF]. Typical functional traits leading to decreased palatability are the presence of thorns (Plate 5) or leathery leaves. Variations in a given interaction between two species can thus arise as a result of changes in the nature and strength of stressors, and the combination of their traits.

Because the dominance of positive or negative interactions in a given plant community may determine strongly its performance and dynamics [START_REF] Callaway | Positive Interactions and Interdependence in Plant Communities[END_REF], many studies have aimed at evaluating the relative importance of these interactions in natural communities [START_REF] Brooker | Facilitation in plant communities: The past, the present, and the future[END_REF]. However, many only consider the variations of pairwise interactions among a limited set of species (e. g., [START_REF] Holthuijzen | Grazing effects on precipitation-driven associations between sagebrush and perennial grasses[END_REF][START_REF] Graff | Testing the role of biotic stress in the stress gradient hypothesis. Processes and patterns in arid rangelands[END_REF], because the number of pairwise experiments required to measure interactions among a diverse set of species is very high. Yet even when they are carried out, pairwise experiments may be insufficient to characterize accurately interactions that occur in a natural community [START_REF] Engel | Can community composition be predicted from pairwise species interactions[END_REF], because of e. g. indirect effects (an interaction between a pair of species is altered by the presence of another). Predictions from the SGH at the community level have also not always matched what is observed at the species-pair level (Soliveres et al., 2014), because a given stressor may have different effects on each species in a community [START_REF] Soliveres | Intransitive competition is widespread in plant communities and maintains their species richness[END_REF]. As a result, to date, no consensus exists on the main drivers of positive and negative interactions in vegetation communities, despite its central role in community assembly and successional dynamics [START_REF] Cramer | What's new about old fields? Land abandonment and ecosystem assembly[END_REF][START_REF] Verdú | Phylogenetic signatures of facilitation and competition in successional communities[END_REF].

Recent work has suggested that the observed co-ocurrence of a pair of species in a given community could be used as a proxy for the effect they have on each other, i. e. their interaction [START_REF] Saiz | Structure and spatial self-organization of semi-arid communities through plant-plant co-occurrence networks[END_REF]. The underlying assumption is that, given a constant abiotic environment, a positive interaction between two species should produce a spatial pattern of aggregation between the individuals of these two species. Conversely, a negative interaction is expected to produce a pattern of segregation. Although co-occurrence-based analyses have shown mixed success in accurately predicting pairwise interactions in general [START_REF] Barner | Fundamental contradictions among observational and experimental estimates of non-trophic species interactions[END_REF][START_REF] Freilich | Species co-occurrence networks: Can they reveal trophic and non-trophic interactions in ecological communities[END_REF], they have been shown

to have some success in predicting community-level properties such as the balance between negative and positive interactions [START_REF] Barner | Fundamental contradictions among observational and experimental estimates of non-trophic species interactions[END_REF]. Some studies have also suggested that associations could highlight shifts in the processes driving community-assembly [START_REF] Lyons | Holocene shifts in the assembly of plant and animal communities implicate human impacts[END_REF][START_REF] Levy | Metabolic modeling of species interaction in the human microbiome elucidates community-level assembly rules[END_REF]. In addition, the equivalence between spatial aggregation and facilitation has been shown experimentally in plant communities [START_REF] Graff | Testing the role of biotic stress in the stress gradient hypothesis. Processes and patterns in arid rangelands[END_REF]. Spatial pattern analyses thus appear as a promising avenue to contribute to the debate over the relative importance intensity of facilitation vs. competition in natural plant communities.

Previous work based on this approach has shown that associations varied according to changes in grazing pressure. In grazed systems, the balance between positive and negative associations can be driven by associative protection, in which a plant's mortality due to grazing is reduced because of unpalatable neighbors. This is expected to produce a unimodal pattern of variation of positive associations [START_REF] Graff | Testing the role of biotic stress in the stress gradient hypothesis. Processes and patterns in arid rangelands[END_REF], where associative protection dominates at intermediate levels of grazing, but is absent when grazing is either low (competi-tion between plants dominate) or very high (disturbance is too high for interactions to be observed, or associative protection disappears).

This pattern has been observed in positive associations [START_REF] Saiz | Changes in semi-arid plant species associations along a livestock grazing gradient[END_REF], but given the fact that associations also tend to vary with aridity [START_REF] Holthuijzen | Grazing effects on precipitation-driven associations between sagebrush and perennial grasses[END_REF][START_REF] Alados | Plant-plant interactions scale up to produce vegetation spatial patterns: The influence of long-and short-term process[END_REF][START_REF] López | Co-occurrence Patterns along a Regional Aridity Gradient of the Subtropical Andes Do Not Support Stress Gradient Hypotheses[END_REF], its generality may be questioned. Other work has also shown that interactions could switch from facilitative to competitive because of grazing, as a result of morphological adaptation [START_REF] Suzuki | Morphological adaptation of a palatable plant to long-term grazing can shift interactions with an unpalatable plant from facilitative to competitive[END_REF], diminishing the potential role of positive interactions in grazed contexts.

No consensus exists yet on how the relative importance of competition and facilitation responds to grazing, leaving a general uncertainty on the effect of stress on the organization of vegetation communities, beyond the sole change in species composition. In addition, despite authors suggesting that facilitation could also be strong in more mesic environments [START_REF] Holmgren | Strong facilitation in mild environments: The stress gradient hypothesis revisited: Strong facilitation in mild environments[END_REF], few studies have focused on the latter. Yet such studies could help understand better the prevalence of facilitation in communities, moving away from nurse-plant asymmetric facilitation and more towards symmetric facilitation (where two plants facilitate each other but have similar lifeforms, [START_REF] Gallien | Frequency and intensity of facilitation reveal opposing patterns along a stress gradient[END_REF], and set a reference for studies carried out in stressed environments.

To contribute to filling this knowledge gap, we investigated the variations of the positive/negative association balance, this time in a more mesic, Mediterranean environment: La Crau, France (mean annual rainfall: 540mm). Contrasting with many previous studies, our work focuses on the variations of associations in herbaceous communities composed in majority of annuals (which are often ignored in association studies, e. g. Soliveres et al., 2014). We investigated whether the variations of positive and negative associations matched previous expectations based on the importance of associative protection, i. e.

whether there was a unimodal relationship in positive associations along the grazing gradient, and how the community-scale balance between positive and negative associations varied along the gradient. We also investigated the species-level association balances (i. e.

whether species tend to engage in more negative or positive associations) to interpret the variations in community-scale trends.

methods description of la crau

La Crau is a Mediterranean semi-arid grassland located in the south of France (Bouches-du-Rhône), whose vegetation has been shaped by sheep grazing for at least 2000 years [START_REF] Henry | The age and history of the French Mediterranean steppe revisited by soil wood charcoal analysis[END_REF]. For this reason, a diverse community of herbaceous plants is maintained (on average 30±12 species/m 2 in this study, mean±s.d.), with a species composition that varies depending on the grazing pressure. All vegetation is herbaceaous with very few perennial higher shrubs, leading to a typical steppe-like visual aspect (Plate 1). Computing species spatial associations requires measuring the spatial structure of plant communities at different levels of grazing pressure. In La Crau, the local grazing pressure (i. e. the number of sheepstime/ha) is known to decrease with the distance from the sheepfold, although its exact value has not been determined for most sites. Previous studies [START_REF] Dureau | Conduite pastorale et répartition de l'avifaune nicheuse des coussouls[END_REF] suggest it decreases as the inverse of distance. In the current absence of a quantitative measure of the grazing pressure, we used an indicator-species-based approach to choose the distances at which the spatial structure was documented. Because previous work identified four vegetation "belts" that differ in species composition around sheepfolds [START_REF] Molinier | La végétation de la Crau[END_REF], we based the spacing of our transects on this classification to make sure we documented a wide range of grazing levels. In a 5x3m area around the quadrats described above, we considered the presence and absence of 16 indicator species known to reflect each vegetation type (taken from [START_REF] Molinier | La végétation de la Crau[END_REF] and [START_REF] Gomila | Relations sol-végétation dans la plain de la Crau[END_REF]). This allowed us to identify in space where each community type occurred. We then laid out two replicate transects in each vegetation type, separated from at least 10m but no more than 75m so that the local vegetation was under similar levels of grazing pressure.

Grazing in La

We documented the spatial structure of vegetation communities using 5 meter-long transects. Each start point and direction of transects were chosen randomly by throwing a pen. For each transect, we laid out a measuring tape, and recorded the length over which every plant individual or part of an individual overlapped to a precision of 2mm (Plate 9 and 10). Positions were determined with an accuracy of 2mm, as estimated from repeated measurements, mainly due to parallax errors when recording plant intersections. For each plant part surveyed, a class of height was also recorded, on a scale ranging from 1 to 5 (classes bounds of 0/2/5/10/20/>20 cm). We also recorded along transects the intersections of pebbles, which drive to some extent the distribution of some species. For example, tillers of Brachypodium retusum tend to resprout from under pebbles, leading to a possible increased spatial association between them. However, we did not include the pebbles in the following analyses as we were interested in the associations between plant species.

We replicated this survey in the surroundings of 6 different sites

(each site comprised one sheepfold). Because in some sites, some veg-etation types ("belts") were absent or impossible to identify, not all belts could be documented for every site. Surveys were carried out during the plant growing season (Apr 15 -June 10) which is the time of the year where 80% of species can be detectable and identified [START_REF] Bourrely | Dynamique actuelle et production primaire nette de l'écosystème steppique de Crau[END_REF], during years 2016 and 2017. This period was 5 included within the grazing season (February -June). Because the structure of plant communities can change along the grazing season, we made sure to distribute the surveys randomly in time.

1 transect pair Figure 3.1: Workflow of the association analysis for community metrics. For each pair of transects, outlying associations are detected, which produces an association matrix of S(S + 1)/2 coefficients for S species, and three raw community-scale metrics K + (percentage of positive links), K -(negative links) and R (association ratio, see main text). The same analysis is carried out for 499 null transect pairs obtained through the randomization of the original transect pair. This creates a null distribution of community-level metrics, which are used to obtain the deviation of observed metrics to what is expected given the same species identities and covers but random spatial structure.

Obs. association matrix

measuring associations

We analyzed species associations using the transect data and a randomization-based approach. This ap-10 proach involves computing a pairwise metric between species on observed data, and comparing it to a distribution obtained from the randomization of the spatial structure.

We computed associations by combining the two replicate transects within each vegetation belt, resulting in one set of associations for each pair of replicate transect. We used the total overlap between all individuals of a pair of species (or between individuals of the same species) along the two replicate transects as a measure of pairwise association. We then randomized the spatial structure by assigning to each individual a new, random position along one of the two replicate transects, and recomputed the overlap metric. By repeating this process 1999 times, this produced, for each species pair, a distribution of overlap metrics that depended only on the abundance of the two species in the dataset (Fig. 3.1). The observed total overlap can then be compared to this distribution.

Associations were classified into positive, negative or neutral by using a cutoff value α applied to the distribution. More specifically, a positive association was considered to occur within a pair of species when 1α/2% of the null distribution was below the observed value. Conversely, a negative association was retained when α/2% of the null distribution was above the observed value (Fig. 3.1). For example, using a classic threshold of 5% (α = 0.05), a positive association would be retained if the observed overlap was above 97.5% of the null distribution. In this study, we chose to display trends based on α = 0.25, but as the choice of this cutoff may affect the results, we carried out a sensitivity analysis to make sure results were robust to a range of thresholds chosen. Given the abundance of all species found in the community and a chosen cutoff value α, this provides information about species that are significantly more frequently found close to each other than expected by chance (positive associations) and species that are less frequently found close to each other than expected by chance (negative associations).

community-scale metrics

For each pair of transects, this method yields a number of positive and negative associations, which can be viewed as a network where nodes constitute species and positive and negative associations are signed links between them. We summarized these networks into three community-scale metrics commonly used in association studies: (i) K + , the percentage of realized positive associations (out of all the S * (S + 1)/2 possible associations that may be realized for S species 1 ), (ii) K -the percentage of realized negative associations and (iii) the relative importance of each association type R = (K + -K -)/(K + + K -) (Saiz et al., 2017; Alados et al., 2017; Ar- royo et al., 2015). R is zero when no positive nor negative associations are detected, and is equal to +1 or -1 when interactions are all positive or negative. In the literature, R has a similar form to the relative interaction index (RII), which has been shown to be a mathematically correct index to measure neighbor-effects [START_REF] Díaz-Sierra | A new family of standardized and symmetric indices for measuring the intensity and importance of plant neighbour effects[END_REF].

The probability of retaining a positive vs. a negative association for a pair of species varies with the total cover of the two species in the transect (Box 3.1), affecting the values of K + , K -and R. The variations of the raw values of these indexes along a gradient may thus reflect variations in species abundances rather than true variations in plant associations. Abundances-corrected values of K + , K -and R need thus to be obtained to correctly interpret the variations of associations along a gradient [START_REF] Pellissier | Comparing species interaction networks along environmental gradients: Networks along environmental gradients[END_REF][START_REF] Tylianakis | Ecological networks across environmental gradients[END_REF].

To do so, we devised a procedure to obtain a distribution of networklevel metrics that reflected expected values given the same identity and species total covers, but random spatial structure. For each pair of transect, we created 199 'null' transects, with randomized positions of individuals along the transects (allowing mixing individuals between replicate transects). These transect pairs have the same species composition and covers, but associations derived from them only reflect the distribution of species covers. We then computed the standardized effect sizes of K + , K -and R compared to this network-level null distribution, K + SES , K - SES and R SES as follows (i.e. Z-score; example given for the latter) :

R SES = R -µ R null σ R null (3.1)
where µ is the mean of the null distribution and σ is its standard deviation.

This allowed us to build a null distribution of community-scale metrics, and disentangle their variations due to changes in species covers from the true effect of the gradient. In what follows, we only present the raw variations of R along the grazing gradient for the Note that this is not S(S -1)/2, because our method also allows the detection of intra-specific associations sake of conciseness, but show the variations of K + SES , K - SES and R SES .

In addition to the variations along the grazing gradient, we investigated the tendency of each species to engage in negative and positive associations. In a similar way to how we derived community-scale metrics, we used, for each species, the percentage of realized positive interactions S + (number of positive associations / total possible number of associations for this species) and negative interactions S -. We then computed the ratio between the two as (S + -S -)/(S + + S -).

As for community metrics, for each species, we corrected these values for species abundances by computing their standardized effect size compared to what was obtained from the null dataset. We compared the SES obtained for each species to the average height of the species (using the midpoints of the height classes) and the average total cover of the species in the transects where it was present.

Box 3.1: Detection bias in plant-plant association analyses

The probability of detecting a significant association (positive or negative) between two species varies with their abundances. To illustrate this, we can consider the simple situation of a transect containing two randomly-placed individuals of different species and identical size s, and consider the association that is retained between the two species depending on the observed overlap. This figure below has been obtained by placing 10,000 time at random two individuals of a given size along a transect, then computing the quantiles of the obtained distribution of overlap. Blue and red bands represent the fraction of the overlap distribution below 2.5% and above the 97.5% quantiles respectively. For a wide range of individual sizes, no negative association can be retained.

In this simple example, negative associations cannot be detected when the cover of the two species is moderately low (below 0.44% of the transect length), regardless of the value of the observed overlap. In other words, even if these two species do not overlap at all, we cannot conclude to a negative association. This arises because the 0.025 (α/2) quantile of the null distribution intersects zero, and the overlap metric cannot be negative (as noted by [START_REF] Saiz | Changes in semi-arid plant species associations along a livestock grazing gradient[END_REF]). As a result, in association networks, all negative associations between rare species are effectively "forbidden links", that can never be derived from data. This has two important consequences on association analyses.

Box 3.1: Detection bias in plant-plant association analyses (cont.)

First, the expected association ratio (R) under random assembly is not neutral (zero) but in fact biased towards positive associations (Figure 3.3), because most natural communities have few abundant species and many rare ones [START_REF] Harte | Maximum Entropy and Ecology: A Theory of Abundance, Distribution, and Energetics. Oxford series in ecology and evolution[END_REF][START_REF] Bazzaz | Plant Species Diversity in Old-Field Successional Ecosystems in Southern Illinois[END_REF]. This could explain the oft-observed dominance of positive associations over negative ones in studies that do not correct for abundance (e. g. Second, because the distribution of abundances of species varies along gradients, the null expectation of associations numbers and ratios will vary accordingly, even in the absence of any biological effect. For example, an increase in rare species along a gradient can shift the association balance towards positive values, even with random spatial structure. It is thus necessary to compare the deviation of these aggregated metrics to what is expected given the distribution of species abundances, rather than their raw values [START_REF] Pellissier | Comparing species interaction networks along environmental gradients: Networks along environmental gradients[END_REF].

regression and sensitivity analyses

We used mixed-effect linear modeling to measure the variations of K + , K -and R (and their abundance-corrected counterparts) along the grazing gradient. We modeled the response of one of the above variables as a linear function of the log of the distance with a normal error structure (e. g. for R, R ∼ a + b * log(dist f ) + ǫ), and considered site-level random effects on the intercept and slope of the regression. We considered an uncorrelated random effect on the slope and the intercept because

given the low number of data points, we could not obtain a relevant estimate of their correlation this correlation was sometimes be 1 and regression diagnostics could not be computed.

We carried out this regression procedure for different values of the association cutoff α and used bootstrapping to derive confidence intervals on the estimate of the slope.

All spatial analyses were carried out in R using custom procedures and the package 'lme4' [START_REF] Bates | Fitting Linear Mixed-Effects Models Using lme4[END_REF] for regressions.

results

Raw variations of R along the grazing gradient were mostly flat and close to zero, suggesting an apparent balance between negative and positive associations (Figure 3.3). This absence of trend was robust to the choice of interaction cutoff α (Figure 3.6) as the effect of distance on the raw ratio was non-significant for a majority of α values. However, null distribution of R values were generally positive along the gradient in all sites, and observed R values almost always below the expected null trends (Figure 3.3, ribbons). In other words, there was always relatively less positive associations than expected given the species richness and covers in transects.

This deviation from random expectation, measured by R SES , changed with decreasing grazing intensity (Figure 3.4). Close to the sheepfold, R SES values were close to the null expectation, while negative values were observed further from the sheepfold. Mixed-effect modelling showed that this decrease was significant for most of the cutoff values α selected (Figure 3.6).

This decrease in R with distance was underpinned by a decrease in the fraction of realized positive associations K + SES , combined with an increase in negative associations K - SES (Figure 3.5). An excess of negative associations was always found, regardless of the grazing level (K - SES values were always positive). Positive associations were close to random in highly-grazed zones, and a deficit was observed in lessgrazed zones (more negative values of K - SES ). These trends were robust for most choices of α cutoff, except for extreme values of α (< 0.10 and > 0.875).

Per-species analyses of associations revealed that most species had a tendency to produce negative associations, as manifested by a negative SES in the association ratio (Figure 3.7). The association ratio of a given species was not significantly correlated to its average height (Kendall's τ = 0.041, P = 0.5) but was negatively correlated with its cover when present in a transect (Kendall's τ = -0.44,P < 0.001) A value close to zero indicates that the distance (i. e., grazing) has no effect on the pattern of association. Negative (resp. positive) values indicate that grazing increases the focal metric (resp. decreases). For example, the raw association ratio R does not significantly vary with grazing for almost all cutoff values (α), but its abundance-corrected version always increases with grazing.

discussion association trends

Using continuous transects documenting the spatial associations of individuals, we found non-random trends in spatial associations along a grazing gradient. Although an excess of negative associations was observed at all grazing levels, very highly grazed areas showed less associations between plants. In area with less grazing, the association balance was highly negative, due to a combined decrease in positive associations and an increase in negative associations.

The classical interpretation of plant-plant associations is that positive associations reflect facilitative processes between plants, while negative associations reflect competitive processes. This hypothesis is consistent with our results, as grazing is known to produce a shift from competitive plant strategies to ruderal plant strategies (le Priol, 2017). Far from the sheepfold, competition processes dominate the assembly of the communities, resulting in an increased number of negative associations. Close to the sheepfold, very high stress due to grazing favors ruderal plants with short lifespans (Chenopodium spp.)

or adapted to frequent disturbance (e. g. Malva sylvestris with a deep taproot, pers. obs), resulting in a spatial structure that reflects less the effect of interactions and more the disruption of spatial patterns by recurrent disturbance.

How do these trends compare to other association-based studies ? We found no unimodal relationship of positive associations along the grazing gradient (contrasting with e. g. [START_REF] Saiz | Changes in semi-arid plant species associations along a livestock grazing gradient[END_REF], although the difference may be explained by methodological differences). This is probably due to the fact that our study is based on the extremely grazed side of the gradient, thus only documenting the collapse of positive and negative interactions under high stress [START_REF] Graff | Changes in sex ratios of a dioecious grass with grazing intensity: The interplay between gender traits, neighbour interactions and spatial patterns[END_REF]. Evidence of associative protection can be observed in the field (Plate 5), but it did not occur frequently enough to be sampled by our systematic survey. The seemingly low presence of facilitative interactions could also be explained by the identity and morphological adaptation of species [START_REF] Suzuki | Morphological adaptation of a palatable plant to long-term grazing can shift interactions with an unpalatable plant from facilitative to competitive[END_REF] We found that the most competitive plants in term of associations were also among the smallest (e. g. Trifolium subterraneum, Figure 3.7),

and there was no significant relationship between the average type of association done by a species and its height. This suggests that the resource under competition in the less-grazed communities is not light (which would be manifested in a tendency for competitive plants to be higher), but rather space and/or local nutrients This interpretation is reinforced by the fact that species with very low association ratios included cespitose grasses (e. g. Dactilys glomerata, Brachypodium retusum), or ground-covering species (e. g. Trifolium subterraneum, Medicago rigidula). Future analyses crossing trait data with associationbased observations would certainly provide insight into the underlying mechanisms driving associations. In addition, our study only looks at above-ground processes, but underground processes could also explain the patterns.

methodological aspects

In this chapter, we showed that strong variations of the association balance along the gradient were masked by variations of species abundances (Figure 3.3). This highlights the need to have community-scale null expectations when interpreting patterns of associations [START_REF] Soliveres | Intransitive competition is widespread in plant communities and maintains their species richness[END_REF]. In fact, since significance of patterns is only tested for the community-scale metrics (R, K + , K -), any α cutoff could be considered at the species level, even one i. e. all associations are included to compute communityscale metrics).

Interestingly, the trends along the grazing gradient varied slightly with the choice of cutoff (Figure 3.6). This calls for systematic testing of the full range of α values (0-1) instead of focusing on a specific values (typically 0.05). In fact, different interaction types may lead to different strengths of associations. For example, using a modelling approach, it can be shown that competition for space produces an excess of weak negative associations (Box 3.2). Using different values as cutoff (or different ranges) may select different types of underlying interaction, whose trends along gradients could be different.

This work shows that association-based approaches can be useful to understand the dominant processes occurring in a community (e. g.

competition, facilitation

). Yet the link between pairwise interactions and spatial patterns remains not fully understood, i. e. which kind of interactions give rise to which kind of spatial pattern. However, the recent rise of interest around co-occurrence methods [START_REF] Morueta-Holme | A network approach for inferring species associations from co-occurrence data[END_REF] and network analyses [START_REF] Graham | Towards a predictive model of species interaction beta diversity[END_REF][START_REF] Delmas | Analysing ecological networks of species interactions[END_REF], along with increasingly-exhaustive trait databases (Tavşano glu & Pausas, 2018) may greatly enhance their potential.

Box 3.2: What interactions are captured by associations ? A modelling approach.

A quick look at the general distribution of association strengths in empirical data reveals that it follows no well-known distribution but has a complex shape. A notable feature of this empirical distribution is the excess of weakly negative association strengths. These associations could possibly represent a widspread and common type of underlying interactions. We used a modelling approach to investigate a possible mechanisms underpinning this pattern.

We devised an individual-based model, in which plant individuals are represented as discs that interact on a 2D plane. The details of this model are presented in Appendix 3. Each individual can affect the mortality and the probability of recruitment of a new individual in its neighborhood. In this model, we also included spatial preemption in part of the simulations. When it was present, propagules cannot access the ground under the canopy of other individuals, so no recruitment of new individuals is possible there. We ran a set of 100 simulations, with and without a space preemption effect, in communities of 20 species until equilibrium. At the end of the simulations, this model produces individuals distributed out in space.

Box 3.2: What interactions are captured by associations ? A modelling approach (cont.). x/y axes represent dimensions in space. A "virtual transect" of length 100 is being represented with a black line.

Because this model is based on disc-like individuals in a continuous space, we measured associations between species by carrying out "virtual transects", which record the range over which individuals intersect a line. This provided us with the same type of data than what was recorded in the field, so we could carry out the same association-based analysis as what was done on empirical data. From these associations, we can then derive the general distributions of associations with, and without space preemption.

Box 3.2: What interactions are captured by associations ? A modelling approach (cont.). Adding the effect of space preemption produces a distribution of association strengths that is strinkingly similar to what is observed in empirical data (compare with Figure 3.2). This suggests that competition for space appears -or another process limiting recruitment or survival under the plant canopies -appears as one of the main process explaining the distribution of observed associations.

T H E R E S I L I E N C E O F I N T E R A C T I O N N E T W O R K S : I N S I G H T F R O M P L A N T -P L A N T A S S O C I AT I O N S context

As shown in Chapter 3, interactions can vary along gradients of biotic stress. In particular, we showed that strong disturbances produced a decrease in competition for space. This chapter focuses on the longterm temporal dynamics of interactions after the release of the biotic stress, i. e. the resilience of interactions. Very few studies have described the temporal variations of interactions (at local scale), and if so it has often only be done for plant-pollinator networks [START_REF] Caradonna | Interaction rewiring and the rapid turnover of plantpollinator networks[END_REF][START_REF] Olesen | Temporal dynamics in a pollination network[END_REF]. This work contributes to filling this gap by describing the temporal dynamics of plant-plant association networks. We focus in particular on a specific structural pattern of association networks: the positive/negative association ratio.

To do so, we adopt a chronosequence approach and document sites that have been previously abandoned at known dates in time. La Crau has a long history of grazing, which allowed us to survey vegetation communities up to approximately 1,600 years after abandonment.

Current sites show a strong, established gradient of grazing centered on the sheepfold. After abandonment, the abandoned sites were still grazed but at much lower levels than originally, and without grazing intensity being spatially centered around the sheepfold. This leaves place for a recovery process of vegetation communities that we document here. This process is remarkably slow: [START_REF] Henry | Origine et dynamique à long terme d'un écosystème herbacé pseudo steppique, le cas de la plain de La Crau[END_REF] main variations of plant-plant associations and frame the chapter in the relevant literature context.

introduction

Ecological systems are increasingly affected by anthropogenic impacts, leaving many communities disturbed, possibly with very different functioning than their pre-disturbance state. This is in particular the case for an increasing number of vegetation communities located in abandoned agricultural fields [START_REF] Cramer | What's new about old fields? Land abandonment and ecosystem assembly[END_REF]. Predicting the future trajectory of this type of communities requires understanding the relative importance of multiple community assembly processes [START_REF] Meiners | Is successional research nearing its climax? New approaches for understanding dynamic communities[END_REF]. These processes determine whether the community will be naturally resilient and recover towards its original, pre-perturbation state [START_REF] Ruprecht | Successfully recovered grassland: A promising example from Romanian old-fields[END_REF], or instead diverge and persist in another state through time [START_REF] Suding | Alternative states and positive feedbacks in restoration ecology[END_REF]. For example, old crops can show an absence of recovery following abandonment, due to competitors impairing the recruitment of the original vegetation [START_REF] Firn | Alternative states models provide an effective framework for invasive species control and restoration of native communities[END_REF]. Disentangling how and where such biotic interactions alter the dynamics of vegetation communities -and their importance relative to other assembly processes -is thus invaluable to anticipate future changes.

Successional dynamics -understood here as the temporal trajectories of communities after a disturbance event [START_REF] Hobbs | Old field succession: Development of concepts[END_REF] -provide an ideal setting to estimate the relative importance of these assembly processes through time after disturbance. A general expectation is that abiotic processes dominate the assembly of the community in the early stages after disturbance, while biotic processes such as competitive exclusion increase in importance in the later stages [START_REF] Connell | Mechanisms of Succession in Natural Communities and Their Role in Community Stability and Organization[END_REF]. However, many exceptions to this general trend can be found (e. g. [START_REF] Verdú | Phylogenetic signatures of facilitation and competition in successional communities[END_REF]. For example, the recovery of many old fields is often limited by the low dispersal ability of native species [START_REF] Fedriani | Assisting seed dispersers to restore oldfields: An individual-based model of the interactions among badgers, foxes and Iberian pear trees[END_REF].

To test the relative importance of different assembly processes, recent correlative approaches have built on functional and phylogenetic community patterns [START_REF] Garnier | Plant functional markers capture ecosystem properties during secondary succession[END_REF][START_REF] Meiners | Is successional research nearing its climax? New approaches for understanding dynamic communities[END_REF]. Where competitive processes are dominant, communities are expected to exhibit overdispersion, i. e. species in the community should be more distinct functionally and phylogenetically than expected through random assembly [START_REF] Webb | Phylogenies and Community Ecology[END_REF][START_REF] Macarthur | The Limiting Similarity, Convergence, and Divergence of Coexisting Species[END_REF]. Functional and phylogenetic overdispersion -hence a possible increased role of inter-specific competition -has been found in late successional stages [START_REF] Purschke | Contrasting changes in taxonomic, phylogenetic and functional diversity during a long-term succession: Insights into assembly processes[END_REF][START_REF] Li | Species colonisation, not competitive exclusion, drives community overdispersion over long-term succession[END_REF]. Yet these approaches have been criticized because the mechanistic bases of overdispersion remain debated [START_REF] Cadotte | Should Environmental Filtering be Abandoned?[END_REF]. These approaches also often ignore the possibility of plant-plant facilitation, which may change the expected successional patterns. For example, [START_REF] Verdú | Phylogenetic signatures of facilitation and competition in successional communities[END_REF] documents a facilitation-based succession which shows phylogenetic overdispersion only in middle successional stages, and not in latestages.

Progress towards a better mechanistic understanding of the drivers of succession -and in particular of the role played by biotic interactions -may come from applying new, independent approaches. In this work, by using a new approach based on plant-plant association networks [START_REF] Saiz | Structure and spatial self-organization of semi-arid communities through plant-plant co-occurrence networks[END_REF], we aim at describing how potential interspecific interactions vary along succession in plant communities.

Association networks describe plant spatial relationships as part of a signed network. When two plants species co-occur more than expected by chance, a positive link is drawn between the two species.

Conversely, a negative link is drawn whenever two species co-occur less than expected by chance. Positive and negative links are thought to reflect on average the sign of the interaction between the pair of species. Thus, the relative proportion of negative and positive links in a given association network could inform on the dominance of positive over negative interactions in a plant community. The variations of this interaction balance during succession may shed a new light on succession mechanisms, and in particular (i) provide an independent alternative to approaches based on functional/phylogenetic relatedness and (ii) include facilitation as well as competition in the relevant processes ocurring along succession [START_REF] Verdú | Phylogenetic signatures of facilitation and competition in successional communities[END_REF].

Association-based approaches are not devoid of criticism [START_REF] Barner | Fundamental contradictions among observational and experimental estimates of non-trophic species interactions[END_REF][START_REF] Freilich | Species co-occurrence networks: Can they reveal trophic and non-trophic interactions in ecological communities[END_REF]. For example, some interactions may produce counter-intuitive associations, such as parasitic interactions, which are negative in essence yet produce a positive association. However, as far as plants are concerned, support for spatial aggregation as the result of facilitation has been shown in experiments [START_REF] Graff | Testing the role of biotic stress in the stress gradient hypothesis. Processes and patterns in arid rangelands[END_REF][START_REF] Valiente-Banuet | Shade as a Cause of the Association Between the Cactus Neobuxbaumia Tetetzo and the Nurse Plant Mimosa Luisana in the Tehuacan Valley, Mexico[END_REF]. Focusing on community-level properties instead of individual pairwise associations may also average out errors at the species-pair level, and provide a more reliable measure of the positive versus negative interaction balance in a given community [START_REF] Saiz | Changes in semi-arid plant species associations along a livestock grazing gradient[END_REF][START_REF] Alados | Plant-plant interactions scale up to produce vegetation spatial patterns: The influence of long-and short-term process[END_REF].

We applied the association network approach to a succession of vegetation communities of the plain of La Crau, a locally semi-arid d). In the early stages after abandonment, the grazing pressure is locally a bit higher than elsewhere as the remains of the sheepfold are still standing and provide shade, where the sheep tend to stay longer (Stage b). We documented the variations of vegetation communities along the line A-B in current and past sheepfolds.

grassland in Mediterranean France. This grassland has been maintained open since the roman era due to a rich sheep husbandry activity [START_REF] Badan | Les bergeries romaines de la Crau d'Arles. Les origines de la transhumance en Provence[END_REF][START_REF] Henry | The age and history of the French Mediterranean steppe revisited by soil wood charcoal analysis[END_REF]. In this plain, current grazing is organized around sheepfolds where sheep stay at night. This results in a local gradient of sheep impact through grazing, nutrient deposition and trampling, whose maximal extension is towards the south-east (the direction under dominant winds). This gradient leads to a local spatial turnover of the species composition. Ruderal plants predominate locally around the sheepfold, within a matrix of a "Coussoul"-type community, dominated by stress-tolerant perennial species such as Brachypodium retusum, Thymus vulgaris and Asphodelus ayardii [START_REF] Molinier | La végétation de la Crau[END_REF]. Through history, some sheepfolds have been abandoned, which results in a locally decreased -but nonzero -grazing pressure and the resorption of more nitrophilous vegetation communities, in favor of the original Coussoul community (Figure 4.3). Since abandonment, a succession-like process has taken place, although it is worth noting that despite abandonment, the vegetation of the sites remains grazed -albeit at much lower intensities.

Because many sheepfolds have a known date of abandonment, the comparison of these sites today provides an ideal setting to adopt a chronosequence approach by space-for-time substitution, allowing us to document the recovery process of the vegetation community through time.

Our study aims at investigating the long-term changes in plantplant associations following grazing abandonment, in an attempt to quantify the resilience of the community, both in terms of species composition and inter-specific interactions.

First, we assessed whether biotic interactions (as measured through associations) varied on the same time scale as the species composition of the community and the abiotic characteristics. To do so, we documented plant associations, species composition and soil attributes of vegetation communities at currently-used and abandoned sites in La Crau, and estimated their trajectory through time by comparing these different sites.

We then focused on testing whether associations supported classic hypotheses regarding succession. Approaches based on functional traits in La Crau (le Priol, 2017) have shown that undisturbed communities tend to contain species with higher leaf dry matter content (LDMC) and lower specific leaf area (SLA), suggesting that the species making up the undisturbed communities may adopt compet-itive strategies [START_REF] Garnier | Inference from Iterative Simulation Using Multiple Sequences[END_REF]. Given these variations, we hypothesized that associations would switch from weak (low number of positive and negative associations) to strongly negative during late successional stages (high number of negative associations).

methods

Surveys and soil sampling

We surveyed 6 currently used sites, 4 sites abandoned in the late XX th century, 3 sites abandoned in the XIX th century and 3 roman sites, abandoned between 88 A.D and 430 A.D (abandonment dates were obtained from [START_REF] Henry | Origine et dynamique à long terme d'un écosystème herbacé pseudo steppique, le cas de la plain de La Crau[END_REF][START_REF] Badan | Les bergeries romaines de la Crau d'Arles. Les origines de la transhumance en Provence[END_REF]. We describe here broadly the survey protocol for the variables used in this chapter, but more details are available in Appendix 2, along with a presentation of the general ecological context of La Crau.

To compute species associations, we laid out transects using a measuring tape, and recorded the beginning and the end of overlap of each plant individual (or part of individual, Plate 9). Transects were laid in pairs -4 pairs in total per site -at different distances from the sheepfold, separated by at least 10m but no more than 75m so that grazing levels were similar. We placed one pair of transect in each one of four community types identified using indicator species (Appendix 2, Figure 5.9). These community types are known to correspond broadly to a certain level of grazing pressure [START_REF] Molinier | La végétation de la Crau[END_REF], so by placing one pair of transect in each type, we made sure that we surveyed a wide range of grazing levels.

We computed the average distance of placement of each transect pair (excluding one site, "Terme Blanc", where the grazing gradient is unusually long), and used this distance to place pairs of transects at abandoned sites.

We merged each pair of transects together to produce one association network per pair. This yielded 4 association networks for each site, at increasing distances from the sheepfold. Because of local disturbances unrelated to grazing, or because indicator species could not be found in the field, some pairs of transects could not be placed, resulting in sites with less than four documented association networks.

For each site, we sampled soil at the average distance of each pair of transect, and complemented this first set of samples with samples taken at quadrat positions. Each soil sample consisted in 200 g taken out of a homogeneous mix of 16 subsamples regularly spaced within a 4 m radius circle (Appendix 2). This yielded between 4 and 6 soil samples per site, which were dried for at least 48h at 40°C, then sent for analyses. Soil analyses included pH (in water), cation total exchange capacity (Metson method), total CaCO 3 , total organic matter, total phosphorus (Olsen method), exchangeable K, Mg, Ca and Na and total Nitrogen.

Selection and derivation of community statistics

We derived for each transect the total cover of each species, then averaged them for each pair of replicate transect. We used this distribution of covers for the pair of transects to compute the number of observed species and the Shannon diversity index. To obtain a univariate value summarizing broad variations in species composition, we first computed a Bray-Curtis dissimilarity matrix between the species compositions of each pair of transects. We then summarized it into a single univariate value by extracting the first axis of a principal coordinate analysis ran on these dissimilarities. The R package 'vegan' v2.5-2 [START_REF] Oksanen | Vegan: Community Ecology Package. R package version 2[END_REF] was used to compute dissimilarities.

We also computed for each pair of transects the average linear plant For the sake of conciseness, we only show the spatial and temporal variations of total nitrogen, which has been shown to be a key factor in driving the plant species composition along succession in La Crau [START_REF] Henry | Origine et dynamique à long terme d'un écosystème herbacé pseudo steppique, le cas de la plain de La Crau[END_REF] and is very strongly correlated to current spatial variations in grazing (Average spearman correlation with distance at current sites: -0.90). We report the variations of other soil variables in

association networks

To obtain the association networks, we computed the total overlap between each species within a pair of transects, then compared this total overlap to what was obtained by randomizing the position of every individual along the transects (and across transects). We retained a positive association when the observed total overlap was within the upper quartile of the null distribution, and a negative one when the total overlap was in the lower quartile (this method is identical to what is described in Chapter 2).

This produced one association network for each pair of transects.

We then summarized these association networks into network-level metrics. First, we computed the percentage of positive and negative associations K + and K -(respectively). We then computed a ratio that quantifies the relative importance of positive over negative associations in a network, R, as (K + -K -)/(K + + K -) (Saiz et al., 2017; Alados et al., 2017; Arroyo et al., 2015). This ration is equal to 1 when all associations are positive, -1 when they are all negative and zero when there is the same number of positive and negative associations in a network.

These three network-level metrics depend on the number of species and their total cover in a given pair of transects. Spatial and temporal variations of these raw numbers may thus reflect these variations in species abundances rather than a true tendency of the community to display a certain type and number of associations (Chapter 3, Box

3.1).

To be able to compare these metrics across gradients, we created 199 'null' transects with a randomized position of each individual (allowing the mixing of individuals across the transects in each pair).

These 'null' transects have a randomized spatial structure, but the same total cover of each species. For each 'null' pair of transects, we computed the network-level metrics, K + , K -and R, which yielded a distribution of 199 null values for each of these metrics. We then computed the standardized effect size of the observed value, compared to its null expectation as (example given for R):

R SES = R -µ R null σ R null (4.1)
where µ is the mean of the null distribution and σ is its standard deviation. We obtained similarly the other network-level metrics K + SES and K - SES . This yields community-level metrics that are broadly independent from the variations of species abundances (Appendix 4, In the following paragraphs, we use y to refer to the property of the community considered, such as species richness, association ratio R SES or the total N measured in a soil sample. We estimated the spatial variations of y at each site by fitting the following saturating function:

y ∼ y 0 * d 1/2 + (y 0 + ∆y) * d d 1/2 + d + ǫ (4.2)
where y is the value of the property considered (e. g. species richness) and d the distance to the sheepfold entrance 1 . In this relation-Why this relationship ? Although the evidence is scarce, grazing pressure P (measured in sheep per unit of area and unit of time) seems to be well approximated by a relationship of the form P = a/(1 + b * d) (Figure 5.5), where a and b are constants and d is the distance to the sheepfold. Assuming the property y of the community under study follows a linear relationship with grazing pressure y = α * P + β, it follows that the relationship between y and d has the form suggested above.

ship, y 0 represents the value of y next to the sheepfold entrance, ∆ y represents the difference between y 0 and the value observed far from the sheepfold (as d → ∞). d 1/2 quantifies the distance at which y = ∆ y /2, i. e. the spatial extent of the variations. ǫ represents a normal error of standard deviation σ.

We used a Bayesian framework to fit this model because it allowed us to include prior information about the general shape of the spatial variations of community properties, and thus maintain parameter estimates within realistic bounds. The details of the model fitting procedure are further described in Appendix 4.

We used 'BayesianTools' v0.1.4 in R (R Core Team, 2018) to fit this model [START_REF] Hartig | BayesianTools: General-Purpose MCMC and SMC Samplers and Tools for Bayesian Statistics[END_REF]. We ran 4 chains for each fit with 100,000 iterations each, using a differential evolution markov chain sampler (DE-MCMC, default sampler in BayesianTools, ter [START_REF] Ter Braak | Differential Evolution Markov Chain with snooker updater and fewer chains[END_REF]. Convergence was assessed visually and by checking that the Gelman-Rubin criterion (Gelman & Rubin, 1992) was under 1.1.

Description of temporal dynamics

To analyze the recovery of community properties over time, we analyzed how the ∆ y values changed with time of abandonment. The value of ∆ y is the amplitude of spatial variations of a given variable around a site: it is expected to be zero when it shows no variation around a site, i. e. a site has fully recovered from its past use for sheep herding. We modelled the variations of ∆ y as a function of time t using the following exponential decay model:

∆ y (t) ∼ ∆ y,∞ + (∆ y,t 0 -∆ y,∞ ) * exp(- log(2)t t 1/2 ) + ǫ (4.3)
In this equation ∆ y,t 0 quantifies the amount of spatial variation observed at current sites, ∆ y,∞ quantifies the long-term spatial variation observed around sheepfolds (it is zero when a site has fully recovered). t 1/2 characterizes the time required for the observed spatial variation ∆ y (t) to be equal to half of its value after infinite time has passed, i. e. (∆ y,t 0 + ∆ y,∞ )/2 -it is thus a measure of the speed of recovery.

We carried out this fit 199 times, using values sampled from the posterior predictive distribution of the previous model. Specifically, we repeated the following process 199 times: for each site, we sam-pled one ∆ y value from the prediction of the spatial variations model, then used these 15 ∆ y values (one per site) to fit the exponential decay model. This allowed us to take into account the uncertainty on spatial variations of community properties into the estimation of temporal dynamics, and estimate confidence intervals around the overall temporal trajectory.

Specific patterns in the variations of associations

Because associations did not show a recovery pattern similar to other aspects of the community, we investigated their variations further. In particular, we investigated two specific patterns: (i) the rapid increase in negative associations close to recently-abandoned sites (abandoned within 9-47.5 years) and (ii) the strong spatial decrease in negative associations in roman sites (1587-1750 years).

We considered the first point by comparing estimates of the percentage of negative associations at the sheepfold (y 0 in equation 4.2)

across classes of time after abandonment (current sites, 0-49.5 years, 119-160 years, 1587-1750 years). We then addressed point two by focusing on the spatial variations of the association ratio, the total nitrogen, species richness, total vegetation cover and relative cover of B.

retusum at roman sites to interpret the long-term trends in plant-plant associations.

Variations of associations with soil nitrogen content

Because patterns suggested a close relationship between nutrient content and the SES of the association ratio, we investigated its general relationship with all sites considered at once. We carried out a piecewise linear regression using one breakpoint to investigate changes in the slope of this relationship. The value of the breakpoint and the coefficients of the linear model were estimated using maximum likelihood and the R package 'segmented' v0.5 [START_REF] Muggeo | Estimating regression models with unknown break-points[END_REF].

results

Spatial variations of community characteristics

The species composition and total nitrogen in the soil showed clear recovery patterns, such that their spatial variations around the sheepfolds decrease with the number of years after abandonment (Figure 4.5) -i. e. over time, assuming a space-for-time substitution.

The temporal trajectories of the amplitude of spatial variations ∆ y confirmed this trend, with species composition and nitrogen showing convergent trajectories of ∆ y towards zero -i. e. over time, the spatial variation of these properties around the sheepfold decreases (Figure 4.4, a-b). This is a very slow process, however, as even after approximately 1600 years after abandonment, a significant spatial variation in species composition can be still observed.

Associations displayed a very different pattern of spatio-temporal variations. The amplitude of variations around a site (∆ y ) did not converge towards zero, but diverged towards a significant non-zero value (Figure 4.4, c). In other words, the spatial variations of associations followed more complex changes over time than a unidirectional recovery pattern. For example, stronger spatial variations were observed around roman sites than in sites of intermediate abandonment The species composition and the total N show a typical recovery pattern, as the spatial variations progressively flatten (red curves) as time from abandonment increases (left to right). On the other hand, the association ratio displays no such recovery pattern, as spatial variations vary and persist through time.

variation of associations at recently-abandoned sites

The spatial variations of associations exhibited rapid changes over time. Most notably, although the association ratio was always negative, its spatial variations around the sheepfold showed a rapid rever- C u rr e n t s it e s 9 4 7 .5 y e a rs 1 1 9 .5 1 6 0 y e a rs 1 5 8 7 1 7 5 0 y e a rs along with the total cover in the transect (-0.11±0.17 % of transect length). This was in favor of an increase of the relative cover of B.

retusum (+0.36±0.23 %), which underpinned an increase in the total aerial volume taken by plants (+1164±277 cm 2 ).

variation of associations with nitrogen levels

The SES of the association ratio was almost always negative regardless of the nitrogen level. However, the piecewise regression identified a turning point (1.83±0.10, mean±s.d., Figure 4.9). For nitrogen levels above this turning point, a decrease in nitrogen yields a decrease in 5 the association ratio (positive slope of 9.8±3.4 SES unit/% Nitrogen).

The reverse was found below this point, as the decrease in nitrogen was associated with an increase in the association ratio value (negative slope of -7.8±3.3 SES unit/% Nitrogen). Total N (g/kg)

Class of abandonment

Figure 4.9: Variations of the SES of the association ratio, as a function of total nitrogen. The line and grey ribbons indicate the results of a piecewise regression (ribbons represent confidence intervals). The black, dashed line has been added by hand and shows a possible expectation of the association ratio in the absence of grazing.

discussion

Observational data from plant communities in La Crau show that the abiotic conditions as well as the species composition followed a recovery trajectory over time after abandonment, although at a very slow pace: 1,600 years after abandonment, these variables are still significantly different near the sheepfold than in their surroundings. This is consistent with previous studies describing the recovery of vegetation communities in La Crau [START_REF] Henry | Origine et dynamique à long terme d'un écosystème herbacé pseudo steppique, le cas de la plain de La Crau[END_REF]. In contrast, plant-plant associations showed a non-recovery pattern. For example, the spatial variations of associations were stronger at the later stages after abandonment (roman sites) than during more recent time periods. This is in stark contrast with our initial hypothesis of a gradual increase of negative associations with time.

It is noteworthy that the spatial variations of associations were similar across sampled sites and time periods. This supports the fact that associations capture a real, information-rich aspect of the organization of the vegetation communities.

Short-term and long-term variations of associations

A first notable pattern captured by associations in the decrease in the association ratio in the early stages after abandonment. This is most likely underpinned by a rebound of competition as grazing is stopped or strongly reduced, while nutrient content remains high. At historical scales however, this is a transient phase and the importance of negative associations stabilizes at a lower level after approx. 100 years (Figure 4.7, g). This rebound of competition reflects fast, possibly annual changes in the organization of vegetation communities.

During spring, at current sites, the impact of sheep around sheepfolds is extremely intense and produces a drastic reduction in cover and vegetation height. However, in the first stages after grazing release, nutrient content is high and a strong and rapid regrowth is often observed of ruderal species. A dominance of negative associations would thus certainly be observed during fall season, resulting in likely strong variations of the ratio during the year (Plate 7).

The plant-plant associations displayed strong variations at roman sites. In particular, these sites had a higher negative association ratio near their entrance than in the nearby "Coussoul" community. The latter was characterized by a dominance of B. retusum, which -de-spite being thought as a competitor -was dominant in a community marked by the absence of negative associations.

This apparent contradiction may be resolved when we consider the nitrogen levels, which were at their lowest in this community. These comparatively low nutrient levels may be what drives the absence of negative associations: in less productive environments, the likely dominant origin of negative associations -competition for space -is lacking and thus the association ratio becomes closer to zero or positive. This depletion in nitrogen could be attributable to the sustained grazing since the roman era, which produces a depletion in nutrients of areas far from sheepfolds, in favor of an enrichment of areas close to them. It could also be a consequence of the larger development of B. retusum, which may preempt nutrients and thus reduce their accessibility to other plants, but this is poorly supported as there is only a weak correlation between the association ratio and the total aerial volume taken by B. retusum (Pearson's r = 0.10, P = 0.58).

Variations of associations and nutrient levels

Be it at recently-abandoned sites or at roman sites, patterns highlighted the close relationship between the association ratio and the nutrient levels available to the community. This relationship showed dual trends. The decrease from very high nitrogen contents to average is associated with an increase in competition -this is because highly nutrient-rich communities also suffer from very high physical disturbance. Were there no grazing but solely nutrient enrichment, the association ratio would have been probably observed as negative (dashed line in Figure 4.9) -reflecting intense competition for space.

However, below a point in nitrogen levels, the association ratio increases towards zero because of a collapse of negative associations, which may be due to the overall decrease in nutrients available to annual plants. It is in these communities that B. retusum increases in biomass, as it may be less affected by the depletion in surface nutrients because of its deeper roots (regularly reaching -40cm). This interpretation remains to be tested, by showing that nutrient limitation is effectively at play in those communities.

Conclusion: the reorganization of communities as identified by associations

To our knowledge, this is the first study applying an association network approach to plant communities outside of arid ecosystems. We show that associations display consistent trends given similar environmental and historical conditions. They may thus be very valuable to measure in situ the changes in the organization of plant communities, and complement previous approaches based on traits and/or phylogenetic relatedness to infer the dominant processes present in a community.

Association analyses confirm the peculiarity of the succession in La Crau following the abandonment of sheepfolds, which converges towards late-successional communities empoverished in species and interactions as nutrient availability decreases. This contrasts with our hypothesis of an increase in competition through time, but it is consistent with the prediction of the stress gradient hypothesis that predicts a decrease in the importance of competitive interactions as productivity decreases. We observed very few positive values of the association ratio, suggesting a relative dominance of competitive processes over facilitation at all stages (but see General Discussion). However, in this study, we used community-level metrics (e. g. the total percentage of negative associations), which may hide variations at the species or species-pair level [START_REF] Sanderson | Patterns in Nature: The Analysis of Species Co-Occurrences[END_REF]. Future developments of this work could extend the analysis down to these levels, along with confirming that nitrogen levels are low enough to effectively produce nutrient-limitation.

C O N C L U S I O N S A N D P E R S P E C T I V E S

Throughout this thesis, we investigated the link between the structure of interaction networks and the response of ecological communities to pertubations. In this section, we summarize and discuss the results of this thesis, then suggest new research directions that have emerged through the development of this work and the many fruitful discussions with colleagues along these past three PhD years.

summary and discussion of research questions

Can active feedback loops be identified from multimodality in diverse ecological communities ? (Chapter 1)

In Chapter 1, we asked whether we could identify vegetation communities in which feedbacks were present, using only a single-time snapshot of their state along a main environmental gradient -water content. To do so, we used an extensive dataset documenting the composition of vegetation communities of subalpine meadows. We found bimodal distributions of species composition in two main areas of the dominant water gradient: dry and wet areas. By reviewing local soil conditions and known previous results about vegetation dynamics in meadows, we concluded that possible active feedbacks could impair the resilience of dry areas in meadows, and provided possible mechanisms. In dry areas, a positive feedback between soil nutrients and vegetation growth may thus exist and be responsible for the bimodality observed. This is a valuable result in a conservation context, where the objective is to identify vegetation communities that may be more vulnerable to disturbance. It is, however, insufficient to characterize precisely the type of dynamics at hand in these fragile communities.

According to the theoretical literature, a strong positive feedback may produce alternative stable states. This could mean that dry communities of subalpine meadows, dominated by Carex filifolia and with coarser soil could be, in fact, stable and persist in time. This stability could be a consequence of erosion or run-off being strong enough in C-filifolia communities to impair the trapping of fine sediments and the regrowth of a more diverse community. However, this cannot be truly demonstrated based on the multi-modality criterion alone, as the presence of multi-modality is not equivalent to the presence of alternative stable states: continuous but sharp responses along an environmental gradient also show increased multi-modality (Box 1.1).

To go further and truly understand the processes behind the bimodal pattern, the most direct way would be to collect temporal data or perform disturbance experiments [START_REF] Petraitis | Multiple Stable States in Natural Ecosystems[END_REF]. However, other correlative approaches may also provide valuable information, such as the analysis of spatial patterns, either through indicators (Chapter

2) or associations (Chapter 3). These spatial patterns help understand the interactions at play in a given community, and could thus also be used to confirm the presence of a feedback in vegetation communities.

How does biotic stress alter plant-plant associations ? (Chapter 3)

We described in Chapter 3 the variations of plant-plant association networks along gradients of biotic stress (grazing). In particular, we focused on an important structural patterns of plant-plant association networks: the number of positive and negative associations, and the ratio between the two. We found a consistent deficit of positive associations, and thus a consistently negative association ratio. This suggests that competition (most likely for space) is the dominant interaction type in these grasslands. This is generally consistent with the fact that the majority of the species surveyed are annuals (55% of the total cover in the dataset), which tend to engage less in specific plant-plant interactions (e. g. allelopathy), other than through the preemption of resources. It is worth noting that this pattern may change in other ecosystems, especially where there may be several different height classes in plants, e. g. where a shrub layer coexist with a lower herbaceous stratum. In previous studies focusing on arid systems [START_REF] Saiz | Redes de Asociación Planta-Planta: Organización de Las Especies y La Comunidad Vegetal, y Adaptación Del Sistema a Factores de Estrés En Ecosistemas Semiáridos[END_REF], Brachypodium retusum engages mostly in positive associations with other plants and higher shrubs in particular. However, it was mostly present in negative associations in the fully-herbaceous communities of La Crau, higlighting a very different behavior of the same species in two different plant communities.

Biotic stress, and more specifically sheep trampling and grazing, reduced the strength and number of associations between species. Interpreting this in terms of interactions, it suggests that under high levels of stress, the dynamics of the plants are governed by the resistance to stress rather than by the interactions they may have with each other. We did not observe any community-level importance of positive interactions such as associative protection, although the latter has been observed in the field (Plate 5). It could play a bigger role in other ecosystems, where differences in palatability between plants are more important, and where grazing levels are not as intense (e. g. [START_REF] Graff | Testing the role of biotic stress in the stress gradient hypothesis. Processes and patterns in arid rangelands[END_REF][START_REF] Saiz | Changes in semi-arid plant species associations along a livestock grazing gradient[END_REF]). This may also be due to the fact that this type of interaction is relatively rare compared to competitive interactions, and thus may be poorly-detectable when considering community-level metrics -such as the association ratio -that average species-specific variations. Extensions of this work would benefit from analyzing these finer-grained variations of associations along the stress gradient, e. g. at the species level or the species-pair level.

More generally, we looked at a single, specific aspect of the association networks we built: the ratio between positive and negative associations. A more thorough analysis of the structure of networks, possibly including the analysis of modularity, nestedness or other structural patterns would probably reveal interesting changes in these patterns.

How resilient are plant-plant association networks ? How do their temporal dynamics compare to other facets of the community ? (Chapter 4)

In Chapter 4, we considered the temporal variations of association networks, using abandoned sites as different snapshots of a given system's recovery process. In terms of network description, we focused on the temporal variations of the balance between negative and positive associations. Interestingly, it proved to be very variable over long time scales and did not show a simple, unidirectional recovery trend such as others properties of the community (e. g. species composition). The resilience of association networks thus does not exhibit a simple, direct return-to-equilibrium response, but rather shows complex variations after a perturbation. These variations were nonetheless consistent across sites of similar historical and environmental condition. This shows that associations describe a real pattern of ecological communities, and contain a lot of information to understand what governs the temporal dynamics of a vegetation community, and in particular the importance of biotic interactions. The association ratio was always mostly negative, highlighting that competition for space was present at every stage of the pseudo-succession. However, variations in this ratio highlighted key phases in the dynamics of the vegetation community -a rebound in competition in the first stages after abandonment, and a decrease of competition in the late successional stages. A general relationship was found between the association ratio and the nutrient levels in the soil. This relationship highlighted a critical nutrient level below which the association ratio becomes closer to zero -i. e. the importance of positive and negative associations is equivalent. At low nutrient levels, competition for space has a reduced importance in structuring the community. This is well in line with predictions from the stress gradient hypothesis (Chapter 3, [START_REF] Bertness | Positive interactions in communities[END_REF][START_REF] Callaway | Positive interactions among alpine plants increase with stress[END_REF], which suggests a decrease of the importance of competition relative to facilitation in less productive environments. This interpretation requires further work to be confirmed. In particular, it may be that these levels of nutrients are low for our dataset, yet not low enough to produce nutrient-limitation in plants.

Plant-plant associations have been mostly used in arid ecosystems.

In these areas, spatial patterns can sometimes be so evident that no statistical analysis is required (e. g. [START_REF] Verdú | The Nested Assembly of Plant Facilitation Networks Prevents Species Extinctions[END_REF], and associations correlate with the already-known importance of facilitation [START_REF] Saiz | Structure and spatial self-organization of semi-arid communities through plant-plant co-occurrence networks[END_REF][START_REF] Motz | Counting Sheep: Roman Settlement and Land Use Patterns in the Crau Plain[END_REF][START_REF] Saiz | The structure of plant spatial association networks is linked to plant diversity in global drylands[END_REF]. For this case, using associations as good proxies for interactions is apparently wellfounded -or at least well-received in the literature.

A major contribution of this PhD is the demostration that associations also show coherent variations in more mesic communities, where facilitation may be less prevalent (Chapter 3 and 4). As such, they appear as a valuable tool to investigate biotic interactions in many vegetation communities. Despite this optimism, it is worth noting that their use is only at its early stages and much work remains to be done to properly understand exactly the information plant-plant spatial associations bring.

perspectives

Towards a better link between interactions and plant-plant associations

The use of plant-plant associations is still in its infancy, as manifested by the low number of publications on plant-plant association networks 1 . A clear upcoming challenge for association-based studies is the lack of a well-established understanding of the link between associations and interactions. While the intuitive equivalence between positive (resp. negative) interactions and spatial aggregation (resp.

segregation) seems well-founded, it still requires a 'leap of faith' at the time of interpreting patterns.

A possibly way forward could be to use independent data such as measurement of traits to make sure the interpretation is correct.

For example, species with competitive strategies could tend to engage more frequently in negative interactions. Some other interactions such as associative protection depend on a specific trait: low palatability, which could correlate with a tendency to produce positive associations. Verifying such hypotheses in the field could guide the interpretation of associations.

Another, more methodological line of research could be to improve the way we compute associations: at this time, they do not take into account many of the complexities found in real-world systems. It is striking that most contributions on plant-plant association networks come from a single author: [START_REF] Saiz | Structure and spatial self-organization of semi-arid communities through plant-plant co-occurrence networks[END_REF][START_REF] Motz | Counting Sheep: Roman Settlement and Land Use Patterns in the Crau Plain[END_REF]; [START_REF] Saiz | The structure of plant spatial association networks is linked to plant diversity in global drylands[END_REF][START_REF] Saiz | Evidence of structural balance in spatial ecological networks[END_REF], and only focus on arid ecosystems

These two examples highlight that interactions may have an effect at different scales and distances to an individual. As a result, point co-occurrence is a too simplistic measure of spatial association, and probably misses interactions between plants. In particular, it may capture better spatial segregation due to competition for space compared to other interactions, which may partially explain the dominance of negative associations in the results presented in this thesis. A possible solution would be to increase the scale of sampling, and consider that co-occurrence happens between two plants when they are closer than a certain target distance (instead of counting co-occurrences only when they overlap). However, this only moves the problem to the choice of this cutoff distance, which is known to affect the results of the analysis [START_REF] Araujo | The geographic scaling of biotic interactions[END_REF][START_REF] Mcnickle | Checkerboard score-area relationships reveal spatial scales of plant community structure[END_REF]. A better method to infer plant-plant interactions should solve this problem and describe aggregation at different scales -e. g. by using the O-ring statistic at different distances. In the next paragaphs, I discuss a possible way to infer better interaction coefficients from the spatial structure of plants.

+ +

a. co-occurrence b. aggregation e.g. shading e.g. associative protection A possible approach could come from individual-based modelling.

We built an individual-based model based on the zone-of-influence principles presented in Figure 5.2 (Appendix 3). This model can produce spatial distributions of plant individuals, given known interaction coefficients between plant species (see also Box 3.2). By 'fitting' this model to data, we could obtain the interaction coefficients (that define curves b, c, d in Figure 5.2) that best reproduce the observed patterns. This is typically done using Approximate Bayesian Computation (ABC) methods [START_REF] Csilléry | Approximate Bayesian Computation (ABC) in practice[END_REF][START_REF] Hartig | Statistical inference for stochastic simulation models -theory and application: Inference for stochastic simulation models[END_REF][START_REF] Hartig | Connecting dynamic vegetation models to data -an inverse perspective: Dynamic vegetation models -an inverse perspective[END_REF].

Such an approach has been applied to derive realistic parameters for individual-based forest models [START_REF] Hartig | Technical Note: Approximate Bayesian parameterization of a process-based tropical forest model[END_REF], but it could be as useful for herbaceous communities. Finally, it is worth keeping in mind that inferring interactions is likely to be a "messy business", involving a synthesis of spatial patterns, trait data and ecological intuition. Bayesian methods have the advantage of explicitely allowing the use of prior, external knowledge, but more likely than not, sound interaction estimations will require combining multiple, correlative and experimental approaches. reorganization, and possibly degradation. In Chapter 2 of this work (Génin et al., 2018a), we focused on indicators of ecosystem degradation. These often rely on the measurement of some kind of aggregation pattern (e. g. spatial autocorrelation), or the size and distribution of vegetation patch sizes, which arise in plant communities because of facilitation [START_REF] Génin | Spatially heterogeneous stressors can alter the performance of indicators of regime shifts[END_REF][START_REF] Schneider | Spatially heterogeneous pressure raises risk of catastrophic shifts[END_REF]. In other words, the spatial association of plants is the basis for indicators of degradation. There is thus a natural connection between these indicators and the use of plant-plant associations. So far, the indicators only consider a binary state (vegetated or not vegetated) and do not use species data. Typically, a single value of vegetation aggregation is derived from the data (e. g. the Moran's I). Association networks could be a way to describe species-specific aggregation in a given landscape. Because pairs of species may react differently to stress, describing the aggregation of species at the species pair level could lead to a more precise and useful indicator of degradation [START_REF] Dakos | Critical slowing down as early warning for the onset of collapse in mutualistic communities[END_REF][START_REF] Holmgren | Strong facilitation in mild environments: The stress gradient hypothesis revisited: Strong facilitation in mild environments[END_REF].

More generally, interaction networks are seldom used as practical management tools, yet they may be very information-rich to identify reorganizations of a community, beyond the sole changes in species composition [START_REF] Memmott | Food webs: A ladder for picking strawberries or a practical tool for practical problems?[END_REF][START_REF] Gray | Ecological networks: The missing links in biomonitoring science[END_REF]. They remain still costly to obtain and complex to analyze, which may be why research has focused more on aggregated properties (biomass, total cover, etc.), yet future research and developments could lower the cost of entry to mapping interactions in empirical communities and increase their use.

The conservation of complex communities

Multiple authors have called for a more holistic vision of what defines an ecological community, that includes not only species but also their interactions (e. g., [START_REF] Tylianakis | Conservation of species interaction networks[END_REF]Memmott et al., 2007a;[START_REF] Gray | Ecological networks: The missing links in biomonitoring science[END_REF] -it is the definition we adopted in this thesis (Figure 0.1).

However, the use of interaction networks in general remains limited:

while common in plant-pollinator studies and food-webs, many other types of interactions lag behind -probably because of the high cost of mapping them. The use of interactions as a unit of description is thus still far from having permeated ecological research.

An interesting lead for interaction research is their use in conservation. In a moving and widely-quoted article, [START_REF] Janzen | The deflowering of central America[END_REF] expressed the problem of the extinction of ecological interactions after deforestation:

« What escapes the eye, however, is a much more insidious kind of extinction: the extinction of ecological interactions. Many of the remaining participants of these interactions will probably hold on for many years, but they constitute little more than a haphazard, semi-self sustaining zoo and botanical garden. Their complex behaviors, morphologies and physiologies seem unnecessarily complex, even maladaptive in the defaunated fencerow-pasture-field habitats they now occupy. »

Decades after this call for the conservation of complex communities rather than just species lists, ecologists are now revisiting central topics and tools of conservation research, but now including interactions as a unit of study; e. g. diversity-area scaling [START_REF] Valiente-Banuet | Beyond species loss: The extinction of ecological interactions in a changing world[END_REF], biotic homogenization (Laliberté & Tylianakis, 2010), beta-diversity [START_REF] Poisot | The dissimilarity of species interaction networks[END_REF]), restoration (Kaiser-Bunbury et al., 2017). However, this mostly concerns plant-pollinator interactions and much work remains to be done for other types of interactions -such as plant-plant interactions. Yet as this work shows, these vary in space and time and some rare and desirable interactions may occur in some places, while the species themselves may not be a conservation target. For example, in La Crau, Brachypodium retusum plays a very original role akin to an "engineering species", yet the species itself is a common and widespread grass, and thus not a conservation

P P E N D I X 2 : D E S C R I P T I O N O F L A C R A U A N D T H E S U R V E Y P R O T O C O L
The plain of the Crau is located in the south of France, and its main area is comprised in a triangle between Arles, Salon and Marseille (Bouches-du-Rhône, 5.3). The Crau area is characterized by a flat landscape where vegetation is dominated by herbaceous species. For this reason, it is often referred to as a steppe (Figure 1), despite not be- In addition, the openness of the area is maintained by a long history of grazing. Dozens of documented roman settlements comprising sheepfolds, along with a rich artifactual and literary evidence supports the presence of a vibrant economy based on sheep herding around La Crau during the roman era [START_REF] Badan | Les bergeries romaines de la Crau d'Arles. Les origines de la transhumance en Provence[END_REF][START_REF] Motz | Counting Sheep: Roman Settlement and Land Use Patterns in the Crau Plain[END_REF].

Charcoal analyses suggest that sheep herding was already active during the Late Neolithic/Bronze Age (approx. -2,000 B.C., [START_REF] Henry | The age and history of the French Mediterranean steppe revisited by soil wood charcoal analysis[END_REF] and remains of enclosures dating from 3,000 B.C. have been found [START_REF] Tatin | Écologie et conservation d'une steppe méditerranéenne la plaine de la Crau[END_REF]. This use of the area for extensive herding persists to this day, but because of historical changes in the sheepherding practices and sheep races, archeological evidence of grazing practices is sparse for some historical periods. Remains of sheepfolds with known abandonment times mostly date from the late roman empire period (100-350 A.D. [START_REF] Badan | Les bergeries romaines de la Crau d'Arles. Les origines de la transhumance en Provence[END_REF], and the XIX th and XX th centurie) [START_REF] Henry | Origine et dynamique à long terme d'un écosystème herbacé pseudo steppique, le cas de la plain de La Crau[END_REF]. Grazing during the middle ages is thought to be the product of a less-specialized economy based on unimproved sheep races that did not require permanent structures [START_REF] Henry | The age and history of the French Mediterranean steppe revisited by soil wood charcoal analysis[END_REF]. In La Crau, flocks of sheep are not confined but driven by shepherds within delimited extensive zones ("Coussouls") of large size (most of them above 70 ha, with a maximum of 540 ha, [START_REF] Tatin | Écologie et conservation d'une steppe méditerranéenne la plaine de la Crau[END_REF].

Most of these zones include a sheepfold built next to a shelter for the shepherd. Sheep are gathered in the sheepfold at night or, because of the large size of the herds (typically 1,600 individuals [START_REF] Tatin | Écologie et conservation d'une steppe méditerranéenne la plaine de la Crau[END_REF], in temporary enclosures next to it. As a result, sheep spend more time near the sheepfold than elsewhere in the area, resulting in a natural gradient of grazing pressure and other impacts related to the presence of sheep (e. g. nutrient enrichement through dung deposition or trampling). This gradient has the greater length towards the south-east, which is the main direction of dominant winds.

Vegetation communities respond to this grazing gradient, resulting in a strong species turnover with ruderal, nitrophilous communities near the sheepfold, and typical "Coussouls"-type communities (dominated by Brachypodium retusum and Asphodelus ayardii) were grazing is at its lowest (Figure 5.9, [START_REF] Molinier | La végétation de la Crau[END_REF][START_REF] Gomila | Relations sol-végétation dans la plain de la Crau[END_REF].

Two other vegetation types (sometimes referred to as "belts"), corresponding to intermediate levels of grazing have been described, one defined by the presence of Onopordum illyricum, and another by the dominance of Trifolium subterraneum. In the absence of accurate measurement of the grazing pressure (but see Figure 5.5), we reused these known vegetation types to guide our surveys and ensure that we documented a large range of grazing pressures. The chapters presented in this thesis are based on a sampling program that spanned three years (2016)(2017)(2018). Surveys were carried out during spring which is the period where sheep are present (typically from March to June [START_REF] Tatin | Écologie et conservation d'une steppe méditerranéenne la plaine de la Crau[END_REF], and when 80% of plants can be identified to the species level [START_REF] Bourrely | Dynamique actuelle et production primaire nette de l'écosystème steppique de Crau[END_REF]. The surveys comprised three different methods in order to document different aspects of plant communities along the direction of maximum extension of the grazing gradient (south-east). This direction had to be adjusted for some sites, where evidence of disturbance was present in the south-east direction (e. g. abandoned fields).

• • • • • • • • • • • • • • • • • • • • • 25
We documented 16 sites spread out in La Crau: 6 which are currently in use, 4 that were abandoned in the second half on the XX t h century, 3 that were abandoned during the first half of the XX t h century or XIX t h century and 3 roman sites (Table 5.1). Abandonment dates were taken from [START_REF] Henry | Origine et dynamique à long terme d'un écosystème herbacé pseudo steppique, le cas de la plain de La Crau[END_REF] and [START_REF] Badan | Les bergeries romaines de la Crau d'Arles. Les origines de la transhumance en Provence[END_REF]. Vegetation communities change in structure and taxonomic composition as spring passes because of the effect of grazing and of the completion of plant life cycles. To make sure this effect was not correlated with any other characteristic of the sites, we randomized during spring the date of survey of each site. We first carried a quadrat-based survey and recorded the presence/absence of all species in a grid of 5 by 5 20x20 cm subquadrats (resulting in a quadrat size of 1m 2 , Figure 5.8). The first quadrat was placed as close as possible to the sheepfold entrance. The spacing between quadrats then followed a geometric series, i. e. the second quadrat was placed at d 2 =2m from the first, then subsequent quadrats were placed such that d n+1 = b * d n were b was between 1.2 and 1.5.

Adjustments to this spacing rule were made to avoid local disturbances unrelated to grazing (Figure 5.7 and 5.9). We also recorded an ocular estimate of vegetation cover and height within each quadrat.

Around each 1m 2 quadrats, we recorded in a 3x6 m area the presence of indicator species (Figure 5.8), which were known to be indicators of the four identified vegetation communities [START_REF] Molinier | La végétation de la Crau[END_REF][START_REF] Gomila | Relations sol-végétation dans la plain de la Crau[END_REF]. These indicator species were then used to ensure that we documented all vegetation types, hence a large range of grazing pressures, and were used to stop the survey once we reached the "Coussouls"-type community.

We recorded the spatial structure of the vegetation, which is required to compute species associations, using 5 meter-long continuous transects. For current sites, we placed a pair of transects in each of the four vegetation types identified by indicator species. Transects were separated by at least 10m but no more than 75 m (to ensure grazing pressure was similar for the two transects, Figure 5.9). This resulted in 8 transects (4 pairs) per site, except when vegetation types could not be observed, resulting in possible missing pairs of transects. We computed the average distance to the sheepfold for each pair of transect, and used this number to place pairs of transects in abandoned sites.

The start and direction of each transect was decided based on a pen throw, avoiding local disturbances such as dirt roads. We then laid out a measuring tape and recorded the length and position of overlap of each part of plant individual on the transect, along with an index of height ranging from one to five (classes bounds of 0/2/5/10/20/>20 cm, Plate 9 and Figure 10). We also recorded the intersection of pebbles along these transects as they can affect the distribution of plant individuals in space.

It is worthy of note that we recorded the change in species composition along gradients twice: with 1m 2 quadrats, then with transects.

In this thesis, we did not use the quadrat data, because transects were Interactions are defined at the species level i. e. all individuals of the same species have the same effect on their neighborhood. An interaction function β i→j (d) describes the effect of species i on species j at a distance d from the center of the individual of species i.

This function is defined as follows: δ here characterizes the distance at which interactions have an effect. The coefficients a ij describe the strength of the interaction between a pair of species and thus define a species-level interaction network. These interactions are directional, i. e. there is no assumption that a ij is equal to a ji , and intra-specific interactions a ii can be included.

β i→j (d) = a ij if d < r i (
In what follows, we use the subscripts r and m on the β i→j and a ij to denote that an interaction affects recruitment (r) or mortality (m).

Temporal dynamics

The model is based on discrete time steps. At time step t, the number and position of newly recruited individuals is determined, as well as the number of dying individuals (removed from the system). The changes are then applied to obtain the system state at time step t + 1.

modelling of plant-plant associations mortality

At each time step, each individual j has a probability of dying P(M j ) (being removed from the system) that is affected by its interactions with all the N other individuals present in the system:

P(M j ) = logit -1   1 + exp β m 0 + N ∑ k=1 k =j β m s k →s j (d jk )    (5.1)
where β m 0 is the log-odds ratio of the base probability p m 0 of dying, log(p m 0 /(1 -p m 0 )). For each propagule spread by existing individuals, a random point

x, y in the system is chosen and the probability P(A x,y ) of the propagule reaching the ground is computed. This probability depends on the k individuals whose canopy intersects the chosen point, and the 'permeability' α s k of those individuals:

P(A x,y ) = ∏ k α s k (5.2)
α s k , the permeability, is a coefficient between zero and one that characterize the probability of a propagule going through the canopy of a given species s k . A higher permeability means that propagule can more easily go through the plant and access the ground. P(A x,y ) effectively quantifies the probability that a propagules goes through all the canopies that intersect the position x, y -and thus access the ground.

If the propagule accesses the ground at location x, y, its probability of growing there P(R j |A x,y ) depends on a base probability p 0 and is modulated by all the N individuals already present in the system at distance d k from the propagule landing point:

P(R j |A x,y ) = logit -1   1 + exp β r 0 + N ∑ k=1 k =j β r s k →s j (d k )    (5.3)
where β r 0 is the log-odds ratio of the base probability p r 0 of recruitment, log(p r 0 /(1 -p r 0 )). To model the fact that plants compete for ressources in addition 5 to space, we introduced a density-dependence on the successful recruitment of new individuals. Each propagule trying to grow in the system has this probability P(R j |A x,y ) further reduced depending of the density d s j (observed cover of species s j ) of the species to which it belongs:

10 P(R ′ j |R j , A x,y ) = P(R j |A x,y )(1 -d s j ) (5.4)
A random uniform number is taken between 0 and 1. If it is below the probability P(R ′ j |R j , A x,y ) then a successful recruitment event takes place (a new individual of species j is added into the system at position x, y).

Simulations

At the start of the simulations, individuals of each species were randomly added in the system until the density of each species reached 1/S (were S is the number of species). At each time step, the variation of species densities were computed by regressing the density of 5 each species against time, using the last 1000 time steps. Simulations were stopped when no species had a temporal variation of its density above 10 -4 .

For the results presented in Box 3.1, we considered two situations.

First, we modelled a situation with no interaction, i. e. all a ij were set 10 to zero. We then considered another situation, but with the physical space occupied by plants taken into account, by setting the permeability of all species α j to zero (where j is the species, Table 5.2). .12: Pair plot of the SES of association ratio with vegetation species richness, total cover of vegetation, bare ground and cover pebble. The total cover of vegetation is given in in % of the length of the transect, and can be above 1 as individuals can overlap.

A P P E N D I X 4 : A D D I T I O N A L M AT E R I A L F O R C H A P T E R 3 A N D 4 5.

correlations between soil variables

Exchangeable Ca (mass %)

rationale for the choice of priors in chapter 4

To estimate the amplitude of spatial variations at each site of community characteristics (e. g. species richness), we fitted a model of the form

y ∼ y 0 * d 1/2 + (y 0 + ∆y) * d d 1/2 + d + ǫ (5.5)
To avoid overfitting and keep the parameters y 0 , d 1/2 and ∆ y , we used a Bayesian framework with weakly informative priors. Because each different community property has different units, the priors could not be set using absolute values, but were set depending on the observed y values at each site.

For each site, sampling yielded a set of values y 1 . . . y n , corresponding to distances d 1 . . . d n to the sheepfold (n is the number of samples taken at a site, typically 2-4 for transects and 4-6 for soil samples).

From this set of points, we extracted y 1 , y n and σ y (the standard deviation of y values). We computed the distance at which the y values are closest to their median x 1/2 , and the average value of y at distances greater than this midpoint y in f . A graphical summary of these variables is present in Figure 5.15.

The priors were then set as follow for the parameters: gie [START_REF] Mouquet | Predictive ecology in a changing world[END_REF].

Si l'on peut raisonnablement identifier les changements d'aire de répartition de chaque espèce à partir de ses besoins abiotiques (e. g. [START_REF] Devictor | Birds are tracking climate warming, but not fast enough[END_REF]), il est difficile de prédire l'effet de perturbations à l'échelle des communautés. En effet, de nombreux systèmes écologiques ne répondent pas linéairement aux pertubations et changements de conditions environmentales, mais peuvent au contraire se réorganiser brutalement, avec parfois des conséquences importantes sur les populations qui dépendent de ces systèmes (Scheffer et al., 2003;[START_REF] Ellis | Long-term effects of a trophic cascade in a large lake ecosystem[END_REF][START_REF] Estes | Trophic Downgrading of Planet Earth[END_REF][START_REF] Vasilakopoulos | Resilience and tipping points of an exploited fish population over six decades[END_REF][START_REF] Lever | The sudden collapse of pollinator communities[END_REF].

Ce type de réponses émerge du fait que les communautés écologiques ne sont pas simples mais complexes, au sens qu'elles sont composées de nombreux éléments en interactions. Les espèces composant une communauté écologique sont liées entre elles par des interactions, qui peuvent parfois être extrèmement nombreuses et diverses. Cette diversité d'interactions rend très difficile la prédiction des effets des perturbations dans les communautés écologiques: par exemple, il est très complexe de prédire les conséquences de l'arrivée d'une nouvelle espèce, comme l'atteste les nombreux "dommages collatéraux" de l'introduction d'espèces dans la lutte contre les ravageurs (e. g., [START_REF] Carvalheiro | Apparent competition can compromise the safety of highly specific biocontrol agents[END_REF].

Malgré cette complexité, les réseaux d'interactions écologiques présentent des patrons structurels qui déterminent la réponse de la communauté entière à des perturbations. Par exemple, certains réseaux d'interactions peuvent présenter des boucles de rétroactions, qui amplifient l'effet des perturbations [START_REF] Wilson | Positive-feedback Switches in Plant Communities[END_REF]. Il a été montré que dans les écosystèmes arides, les interactions facilitatrices entre plantes peuvent conduire à des réorganisations abruptes suite à un changement d'aridité, même faible [START_REF] Kéfi | Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems[END_REF][START_REF] Rietkerk | Self-organized patchiness and catastrophic shifts in ecosystems[END_REF]. Dans ces systèmes, les plantes poussent en îlots (patch), ce qui permet de maintenir une certaine disponibilité en ressources localement. Lorsque le stress augmente, la présence des interactions facilitatrices et la structure en patch permet d'absorber l'augmentation du stress, jusqu'à un certain point critique où la communauté se dégrade et devient désertique.

Les interactions présentes dans les systèmes écologiques ne conduisent pas toujours à de telles réorganisations abruptes. D'autres études montrent également que la structure de nombreux réseaux d'interactions entre espèces aurait tendance à augmenter la stabilité de la communauté [START_REF] Lever | The sudden collapse of pollinator communities[END_REF][START_REF] Verdú | The Nested Assembly of Plant Facilitation Networks Prevents Species Extinctions[END_REF][START_REF] Memmott | Tolerance of pollination networks to species extinctions[END_REF][START_REF] Stouffer | Compartmentalization increases food-web persistence[END_REF][START_REF] Thebault | Stability of Ecological Communities and the Architecture of Mutualistic and Trophic Networks[END_REF].

De nombreux travaux ont donc été dédiés a l'identification de la Les nombreux travaux portant sur les réseaux d'interactions plantespollinisateurs ont montré que leur structure déterminait la réponse de la communauté aux perturbations [START_REF] Memmott | Tolerance of pollination networks to species extinctions[END_REF], mais l'intensité des interactions elles-mêmes était également affectée par les pertubations (Memmott et al., 2007b). Pour décrire la dynamique Cover (%) .

-10 -5 0 . répondre aux perturbations. Une dominance de la facilitation entre plante à l'échelle de la communauté peut potentiellement produire des boucles de rétroactions [START_REF] Kéfi | When can positive interactions cause alternative stable states in ecosystems?[END_REF], ou encore permettre à un nombre plus grand d'espèces de coexister [START_REF] Saiz | The structure of plant spatial association networks is linked to plant diversity in global drylands[END_REF]. abstract Ecological systems are not simple but composed of many different elements (species, for example) interacting with each other. These networks of interactions exhibit structural properties that determine ecological systems' ability to absorb and recover from perturbations. Mapping interactions along with their changes in time and space is therefore key to understand and predict empirical communities' response to global changes. In this thesis, we used plant communities as model systems (i) to explore how spatial patterns may help identify feedbacks loops which make communities more fragile to upcoming changes and (ii) to map species interactions in empirical communities and describe how they change along stress gradients and recover from perturbations. To do so, we used two datasets documenting plant communities in subalpine meadows (USA) and Mediterranean grasslands (France).

Est. di

Les communautés végétales

Our results show that feedback loops can be inferred to some extent from the spatial patterns of plant communities and hence help identify communities that may respond more abruptly to perturbations. Going to a more detailed level of description, plant-plant interactions (as measured through spatial associations) were shown to respond strongly and consistently to stress but exhibited a weak resilience to disturbances. This work shows that plant-plant interactions -which are linked to the response of the community to perturbations -can be uncovered using spatial patterns. It paves the way towards a better understanding and a better anticipation capacity of how ecological communities might reorganize when subject to disturbances.

résumé Les systèmes écologiques sont complexes car composés d'une multitude d'éléments en interaction. Ces interactions, entre espèces par exemple, forment des réseaux qui présentent des propriétés structurelles déterminantes pour la réponse du système écologique entier aux perturbations. Pour mieux identifier cette réponse, il est donc important de cartographier les interactions présentes dans les communautés écologiques et de comprendre leurs variations dans le temps et l'espace. Dans ce travail, nous avons utilisé les communautés de plantes comme systèmes écologiques modèles afin (i) d'identifier à partir de patrons spatiaux certains motifs présents dans les réseaux d'interaction écologiques (les boucles de rétroaction) et (ii) de cartographier les réseaux d'interactions (tels que mesurés par les associations spatiales entre plantes) le long de gradients de stress. Pour ce faire, nous avons utilisé deux jeux de données documentant des communautés de plantes de clairières subalpines (Etats-Unis) et méditerranéeennes (La Crau, France). Nos résultats montrent que les boucles de rétroaction peuvent être inferées à partir des patrons spatiaux présents dans les communautés de plantes, permettant ainsi d'identifier des communautés pouvant répondre de manière abrupte aux perturbations. Les interactions entre plantes (déduites de leurs associations spatiales) dépendent fortement du stress appliqué à la communauté, et présentent une résilience faible aux évenements de perturbation. Ce travail montre que les interactions entre plantes peuvent être cartographiées in situ à partir des associations spatiales. Il ouvre la voie vers une meilleure compréhension et capacité d'anticipation de la réponse des communautés écologiques face aux perturbations.

1. 1

 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 2 fostering the detection of ecosystem shifts by 15 implementing spatial early-warning signals 2.1 Background . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 0 . 4

 04 Figure 0.1 Interplay between interaction network, resilience and disturbance . . . . . . . . . . . . . . . . . . Figure 0.2 Interaction networks in ecological systems . . Figure 0.3The stable equilibria of the overgrazing model 5

Figure 0. 5

 5 Figure 0.5 Types of networks found in co-occurrence analyses . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1. 1

 1 Figure 1.1 Picture of a subalpine meadow . . . . . . . . . Figure 1.2 Picture of a dry and wet meadow plant com-10

Figure 1. 3

 3 Figure 1.3 General view in a meadow . . . . . . . . . . . .

Figure 1. 4

 4 Figure 1.4 Principle of component analysis on synthetic data . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1

 1 Figure 1.5 Detrended Correspondence Analysis results .15

Figure 1. 7

 7 Figure 1.7 Detrended Correspondence Analysis results on the wet meadow subset . . . . . . . . . . . . . .

Figure 1. 8 Figure 2

 82 Figure 1.8 Detrended Correspondence Analysis results on the dry meadow subset . . . . . . . . . . . . . .

Figure 2. 2

 2 Figure 2.2 Example Generic EWS trends . . . . . . . . . .

Figure 2. 3

 3 Figure 2.3Example Spectral Density Ratio trends . . . . .

Figure 2. 4

 4 Figure 2.4 Patch-based indicator trends . . . . . . . . . . .

Figure 2. 5 Figure 3 . 1

 531 Figure 2.5 Individual patch-size distribution fits . . . . . 25

Figure 3. 4 Figure 3

 43 Figure 3.4 SES of the association ratio along the grazing gradient. . . . . . . . . . . . . . . . . . . . . . .

Figure 3. 6

 6 Figure 3.6 Regression and sensitivity analyses of the variation of associations . . . . . . . . . . . . . . .

Figure 3. 7

 7 Figure 3.7 Associations per species . . . . . . . . . . . . .35

Figure 4. 3

 3 Figure 4.3 Response of vegetation communities after sheepfold abandonment . . . . . . . . . . . . . . . . .

Figure 4. 4 Figure 4 . 5

 445 Figure 4.4 Temporal trajectories of community characteristics . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4. 6

 6 Figure 4.6 Estimates of the % of negative associations at the sheepfold (SES) . . . . . . . . . . . . . . . .

Figure 4. 7

 7 Figure 4.7 Spatial variation of positive, negative associa-10

Figure 4. 8

 8 Figure 4.8 Spatial variations of vegetation communities at roman sites . . . . . . . . . . . . . . . . . . . . .

Figure 4. 9

 9 Figure 4.9 Variations of the SES of the association ratio with nutrient levels . . . . . . . . . . . . . . . .

Figure 5. 2

 2 Figure 5.2 The effect of interactions on the probability of occurrence . . . . . . . . . . . . . . . . . . . . .

Figure 5. 3

 3 Figure 5.3 General situation of the Crau. The main area is 20

Figure 5. 4

 4 Figure 5.4 Aerial view of a sheepfold and sheep daytime trajectory . . . . . . . . . . . . . . . . . . . . . .

Figure 5. 5 Figure 5 . 6

 556 Figure 5.5 Grazing pressure as a function of distance to the sheepfold . . . . . . . . . . . . . . . . . . . .

Figure 5. 8 Figure 5 . 9

 859 Figure 5.8 Information recorded for each quadrat survey and details of soil sampling . . . . . . . . . . .

Figure 5. 10

 10 Figure 5.10 Model representation of interactions . . . . . .

Figure 5. 11

 11 Figure 5.11 Life cycle of an individual and interactions . .

Figure 5. 12

 12 Figure 5.12 Correlation of the SES of association ratio with 35

Figure 5. 14

 14 Figure 5.14 Pair plot of relationships and correlation between soil variables . . . . . . . . . . . . . . . .

Figure 5. 15

 15 Figure 5.15 Graphical representation of the parameters used for priors . . . . . . . . . . . . . . . . . . . . . . 152 Figure 5.16 La relation entre réseau d'interaction, résilience et perturbations . . . . . . . . . . . . . . . . . . 155 Figure 5.17 Bimodalité dans les communautés sèches des 5

Figure 0 . 1 :Figure 0 . 2 :

 0102 Figure 0.1: An ecological community is defined by a set of species and their interactions. These interactions can show structural patterns that determine the resilience and stability properties. The latter, combined with the characteristics of disturbance events, affect back the species abundances and interactions, effectively changing the original ecological community.

Figure 0 . 3 :

 03 Figure 0.3: Bifurcation diagram of the overgrazing model (with r=1, K=10, V 0 =1), along a gradient of c values.Continuous lines indicate stable states, and the dashed line indicates the limit of the two basin of attractions where multiple stable states are present for the same grazing values. This line delineates two basins of attraction, which are the set of values from which a system will converge to a particular stable state (grey arrows). White dots indicate critical points, at which one of the multiple stable states becomes unstable.

Figure 0 . 4 :

 04 Figure 0.4: Potential curves of the overgrazing model (red curves). Stable states in the bifurcation diagram correspond to the bottom of valleys in the potential curves (black and grey dots). As the environmental conditions change (here, the grazing c), the shape of the curves change and the position and number of stable states change.

Box 1 :

 1 An example null-model-based analysis of species co-occurrences (cont.) Distribution of the null number of co-ocurrences for a pair of species. If the observed number of co-occurrences is above the 97.5% quantile of the null distribution of co-occurrences (red line), a positive association is retained. A similar test yields negative associations in the lower tail of the null distribution (blue line).

Figure 0 . 5 :

 05 Figure 0.5: Types of networks found in co-occurrence analyses. Networks can represent directional interactions (a,b,c, e. g. food-webs), or undirectional interactions (d,e,f). The same classification applies to association networks derived from co-occurrence analyses, which represent associations as signed, qualitative directed or undirected links (b, e). Weighted networks can be also obtained, in which each association has a given strength (c,f). In most studies, the networks obtained through co-ocurrence analyses are signed and undirected (i. e. of type e or f).

Figure 1 . 1 :

 11 Figure 1.1: Subalpine meadow often appear as clearings within a continuous pine forest. They occur where water accumulates and favors a diverse herbaceaous community. Photo: Steven Lee (USGS)

Figure 1 . 2 :

 12 Figure 1.2: The species composition of meadow vegetation communities responds first and foremost to changes in local water availability. The photo on top pictures a community of Carex in a regularlyflooded meadow area, while the bottom photo shows a community of Carex (a different species) in dry areas of meadows. Photos: Steven Lee (USGS)

  For a given system, identifying changes in the number of modes -the modality (Petraitis, 2013) -of the empirical distribution of state values may thus identify where non-linear responses occur along a gradient. These areas may be of special interest to carry out further research or experiments, as perturbations or small changes in environmental conditions may have there a larger effect than expected (Box 1.1).

Box 1. 1 :

 1 The link between modality and non-linear responses

Figure 1 . 4 :

 14 Figure 1.4: Principle of component analysis on synthetic data with addednoise. On univariate data (a,c), slicing the density yields a single component for a majority of chosen thresholds. On bimodal data, the same procedure yields two components for most thresholds chosen (b, d). Note that very low thresholds identify spurious components because of outliers in the density, but this effect disappears quickly as higher thresholds are used

  Unfortunately not shown in this thesis.

Figure 1 . 5 :Figure 1 . 6 :

 1516 Figure 1.5: Detrended Correspondence Analysis results on the subalpine dataset (a). Species labels indicate where individual species are dominant, with label size proportional to the overall abundance of the species in the dataset. The bottom panels (b1-4) show the same graphs, but only with points comprised within different ranges of water content. Increased modality (i. e. multimodality) can be observed for water contents below -1 s.d. and above +1 s.d.

Figure 1 Figure 1

 11 Figure 1.7: (a) Detrended Correspondence Analysis results on the wet meadow subset. Species labels indicate where individual species are dominant, with label size proportional to the total abundance of the species in the subset. A density estimate is added on top of the points. (b) Species composition and differences in plot attributes between the two groups, stars represent significance of a Mann-Whitney test with P-values corrected for multiple tests using the Benjamini-Hochberg procedure (one star for P < 0.05 and two for P < 0.01). The grey bars show the estimated difference in (pseudo)medians between the two groups.

  community contrasting with a more diverse community of mesic species. The C. vesicaria-dominated community had a larger amount of flooded or saturated soil and more silt, which are typical characteristics of streamside flooded pools with standing water (Baldwin & Hickman, 2012). The mesic community occurring for similar water contents was represented by a more diverse species composition that occurred outside of the areas with standing water. Here, species composition maps well onto independent, external environmental characteristics linked with the absence or presence of standing water. A possible transition from one community type to the other would likely be the result of a major modification of the local water regime.

  (i) under low or no stress, a spanning cluster of occupied(vegetated) cells is present. (ii) As stress increases, this single patch breaks down and the PSD resembles a power-law (x -α ). (iii) As stress increases further, vegetation patches break down and the PSD fits best to a truncated power-law (x -α e -βx ). (iv) In the final stages before reaching a fully empty state, only small patches persist and the patch size distribution is best-described by an exponential distribution (e -βx ) with presence of a spanning cluster of empty cells. Identifying where in these four phases an ecological system lies could thus constitute an indicator of degradation level[START_REF] Kéfi | Robust scaling in ecosystems and the meltdown of patch size distributions before extinction: Patch size distributions towards extinction[END_REF], but see[START_REF] Schneider | Spatially heterogeneous pressure raises risk of catastrophic shifts[END_REF];[START_REF] Sankaran | Patchiness and scale-free correlations: Characterising criticality in ecosystems[END_REF]).

  workflow allows computing EWS from raster data, typically from re-

  10 their significance compared to a random spatial structure, and (iv) display the results (Figure 2.1).

Figure 2

 2 Figure 2.1: 'spatialwarnings' workflow. Black headers indicate the relevant package functions for each step.

(

  2017) on remote sensing imagery. In this study, the authors used remotely-sensed images along gradients of increasing annual rainfall. The aerial images, stored in the R object serengeti were preliminarily classified so that TRUE values (pixels in the image) represent grassland areas and FALSE values represent woodland areas. The transects document the shift of the ecosystem from a grassland to a woodland state around an annual rainfall of 730mm/year. The object serengeti.rain contains the rainfall values for each matrix. Generic EWS -variance, skewness and near-neighbor correlation (as measured by Moran's I) -can be computed using the generic_sews() function. However, the computation of these indicators differs depending on the data type.When working with binary data (TRUE/FALSE values), variance and skewness only reflect the proportion of TRUE values in the data and do not depend on the spatial structure[START_REF] Kéfi | Early Warning Signals of Ecological Transitions: Methods for Spatial Patterns[END_REF]. Coarsegraining has been shown to be a useful procedure to make these indicators actually reflect the spatial patterns(Sankaran et al., 2017a). This procedure divides the input matrix of size N × N (for a square matrix) into submatrices of size s × s. The average of the pixels is taken in each submatrix, thus producing a N/s × N/s matrix with continuous values. Even when a matrix already contains continuous values, in some cases, coarse-graining may nevertheless be required to correctly compute EWS, and is therefore provided optionally.Sankaran et al. (2017a) provides more details on the procedure and on the choice of coarse-graining length.Note that, in the package, coarse-graining is only included by default in the computation of spatial variance and spatial skewness. Nearneighbor correlation can be forced to be computed on coarse-grained data using the moranI_coarse_grain option. Using a subsize s = 5 as per[START_REF] Eby | Alternative stable states and spatial indicators of critical slowing down along a spatial gradient in a savanna ecosystem[END_REF], we can reproduce their analysis using the following code: serengeti.ic <-generic_sews(serengeti, subsize = 5, moranI_coarse_grain = TRUE) summary(serengeti.ic) The significance of generic EWS can be assessed by comparing the observed indicator value to a distribution of values obtained from matrices with a randomized spatial structure. The generic function indictest() can perform this operation by permuting the position of the cells in the matrix a great number of times (999 by default). To carry out this operation efficiently, much of the package is implemented in compiled code, and the package will use the global option mc.cores to carry out parallel processing. serengeti.test <-indictest(serengeti.ic) 5 summary(serengeti.test) Calling the plot() function on serengeti.test displays the results of the computation (Figure 2.2):

MoranFigure 2 . 2 :Figure 2 . 3 :

 2223 Figure2.2: Generic EWS trends for the serengeti dataset[START_REF] Eby | Alternative stable states and spatial indicators of critical slowing down along a spatial gradient in a savanna ecosystem[END_REF] along an annual rainfall gradient. This figure is produced by calling the plot() function on the output of indictest(). First row: mean value of each input matrix (here average grassland cover). 2nd-4th rows: the three generic EWS -variance, skewness and near-neighbor correlation (Moran´s I). Black lines represent the observed values, and gray ribbons the 5%-95% quantiles of the null values (obtained by randomizing the spatial structure). The red line displays the value at which the ecosystem shifts from savannah to forest (added using ggplot2; see Appendix 1). The three indicators are found to increase prior to the shift.

Figure 2 . 4 :

 24 Figure 2.4: Outcome of the plot() function called on the output of patchdistr_sews().The upper panel indicates the best fitted PSD and shows whether there is percolation (presence of a spanning cluster) of full (TRUE) or empty cells. Mean cover is displayed as a continuous line. As stress increases, best-fit of PSD switches from a power-law (pl) with percolation of full cells, followed by a truncated power-law (tpl) and an exponential (exp) with percolation of empty cells. The bottom panel shows PLR decreasing along the stress gradient.

Figure 2 . 5 :

 25 Figure 2.5: Individual PSD fitted to each matrix of the forestgap dataset. As stress increases (top to bottom and left to right), the truncation of the PSD increases and the shape of the best-fitting distribution goes from a power-law (pl), to a truncated power-law (tpl) and an exponential (exp). Along the same gradient, the PLR decreases (represented by the blue double-arrowed bar on top of each panel).

  why, since the late 90's and early 2000's, much work has focused on measuring the relative importance of competition vs facilitation in vegetation communities. We contribute to this debate using a relatively new approach based on co-occurrence networks, and describe the variations of the competition/facilitation balance in vegetation communities of La Crau (France) along grazing gradients. This work is based on vegetation surveys carried yearly during along the course of this PhD program (springs of 2016, 2017 and 2018, and fall of 2017 for soil sampling), totalling an approximate amount of 4 months of field work. It was carried out in collaboration with Thierry Dutoit at IUT d'Avignon, along with Axel Wolff and Laurent Tatin from the Réserve Naturelle des Coussouls de Crau.

  Crau is centered around numerous sheepfolds (pens for sheeps) spread out in the steppe, where sheeps are gathered at night. Sheep activity is most intense around the sheepfold and decreases with the distance from the sheepfold. This results in local changes in the identity and diversity of vegetation communities in an otherwise continuous "Coussoul"-type community[START_REF] Molinier | La végétation de la Crau[END_REF], dominated by the perennial grass Brachypodium retusum (Plate 4). Ruderal nitrophilous plants are most abundant close to sheepfolds, and decrease in abundance as distance to the sheepfold entrance increases. This very-well defined gradient constitutes an ideal setting to test the influence of grazing (and other stresses associated with grazing, such as trampling and redistribution of nutrients) on the associations of vegetation communities.surveysWe selected six sites, i. e. six sheepfolds, for sampling that were spread out in the studied area (Figure5.6), and displayed typical vegetation turnover as a response to grazing.We first carried out a quadrat-based survey to document the variations of species covers along the grazing gradient. Each quadrat had an area of 1x1m and was subdivided into a grid of 25 subquadrats (20x20cm each), in which the presence/absence of each species was recorded. The first quadrat was placed as close as possible to the entrance of the sheepfold, avoiding any local disturbance not related to grazing (e.g. presence of a dirt road). The spacing between quadrats followed a geometric series, i. e. the second quadrat was placed at d 2 =2m from the first in a general south-east direction (maximal direction of extension of the gradient), then subsequent quadrats were placed such that d n+1 = b * d n were b was between 1.2 and 1.5 depending on the site (Appendix 2). Adjustments in the quadrat spacing or in the general direction of the transect were made to avoid local disturbances unrelated to grazing.

Figure 3 . 2

 32 Figure3.2 Type of association deduced from the observed overlap between two individuals of a given size (for α = 0.05). Blue and red bands represent the fraction of the overlap distribution below 2.5% and above the 97.5% quantiles respectively. For a wide range of individual sizes, no negative association can be retained.

Figure 3 . 3 :

 33 Figure 3.3: Variations of the raw association ratios R and its expected valuegiven the species covers along the grazing gradient (measured by the distance to the sheepfold). Each panel corresponds to one surveyed site (along the rows). The raw association ratio showed mostly flat variations along the gradient, with values close to zero (same number of positive vs negative associations). The expected values of the ratio given the species covers and abundances were always positive.

Figure 3 . 4 :

 34 Figure 3.4: Variations in the SES of the association ratio R + SES along the grazing gradient. A value close to zero indicates that the ratio value is fully explained by the species composition and covers of the given pair of transect. Negative values indicate a lower association ratio that what is expected given the species covers.

Figure 3 . 5 :Figure 3 . 6 :

 3536 Figure 3.5: Variations in the SES of percentage of realized positive associations (K + SES ) and negative associations (K - SES ) along the grazing gradient. A value close to zero indicates that the observed number of negative and positive associations is close to its expected values given the species composition and covers. Negative (resp positive) values indicate a deficit (resp an excess) in the given type of interaction.

Figure 3

 3 Figure 3.7: (a) SES of association ratio for each species (only species with at least 200 potential associations were included in this graph (see species codes in Appendix 1). (b) Relationship between the association ratio and the average height of the plant. (c) Relationship between the association ratio and the average cover of the plant in the transects where they are present.

Figure 3 . 8

 38 Figure 3.8 Empirical distribution of 1024 randomly-selected association strengths (i. e. the deviation of the observed total overlap between pairs of species to the mean of the random expectation, divided by the s.d. of the random expectation) The y-scale is non-linear (square-root transformation).

Figure 3 . 9

 39 Figure 3.9 Example output of the model. Discs are individuals of different species (color).x/y axes represent dimensions in space. A "virtual transect" of length 100 is being represented with a black line.

Figure 3 . 10

 310 Figure 3.10 Observed distribution of 4096 association strengths from model data, with and without the effect of space preemption.

Figure 4 . 1 :

 41 Figure 4.1: The sheepfold Guillaume Orcel, which has been abandoned between 1968 and 1971 (Henry, 2009).

Figure 4 . 2 :

 42 Figure 4.2: Artist depiction of a Roman sheepfold (image taken from a interpretative sign in Peau de Meau). See also Plate 8.

Figure 4 . 3 :

 43 Figure 4.3: Response of vegetation communities at a currently-used site, recovery through time, and hypothesized spatial distribution of grazing pressure. A currently-used site will show a strong local variation of grazing pressure (Stage a), which results in a local change of vegetation communities (colors). After abandonment, a "rewinding" process occurs and the "Coussoul"-community (grey) colonizes back the disturbed areas (Stages b-d). In the early stages after abandonment, the grazing pressure is locally a bit higher than elsewhere as the remains of the sheepfold are still standing and provide shade, where the sheep tend to stay longer (Stage b). We documented the variations of vegetation communities along the line A-B in current and past sheepfolds.

  volume, i. e. the sum of the size of each individual (or part of individual, Plate 9) multiplied by its height, to obtain an estimation of the above-ground plant community volume. The volume taken by an individual was computed as follows: if individual i overlaps the transect between a i and b i , and has a height h i , then its volume is (b i -a i ) * h i . Because we only recorded classes of heights (1/2/3/4/5) rather than actual height in the field, we carried out the average using the midpoints of each height class (1/3.5/7.5/16cm and 25cm for class 5).

Appendix 4 .

 4 It is worth noting that most soil variables are correlated with each other and in particular with the total nitrogen content (average Pearson's r=0.67±0.2 s.d.), thus yielding similar trends in space and time (Appendix 4, Figure 5.13 and 5.14).

Figure 5 .

 5 Figure 5.12), and thus allows meaningful comparisons of SES values

(Figure 4 . 4 :Figure 4 . 5 :

 4445 Figure 4.4: Temporal trajectory of the spatial variations of species composition, the soil total Nitrogen and the SES of the association ratio. The trend line corresponds to the exponential decay model described in Equation 4.3, and the ribbon indicate 95% confidence interval around the predicted mean value .

5

  sal after abandonment(Figure 4.7,. This reversal was the consequence of a strong and significant increase in the dominance of negative associations close to the sheepfold (Figure4.7, a-b, and Figure4.6).

Figure 4 . 6 :Figure 4 . 7 :

 4647 Figure 4.6: Estimate of the SES of the % of negative associations at the sheepfold, for each class of abandonment. Letters on top indicate significantly different groups, based on pairwise Wilcoxon tests and a significance threshold of 0.10.

Figure 4 . 8 :

 48 Figure 4.8: Spatial variations, at roman sites, of the association ratio, total nitrogen, species richness, total cover of plants (in % of the transect length), and relative cover of B. retusum.

  For example, co-occurrence can only imperfectly detect aggregation. Here, we considered that two individuals overlapped when they occupied the same point along a linear transect (i. e. carrying out the continuous equivalent of point-intersect sampling). In fact, the zone of influence of many plants extends beyond their exact location in space and some plants can be spatially aggregated yet never totally overlap (Figure5.1). For example, a plant can provide shade and improve local conditions for another one, yet never occur at the same exact position in space[START_REF] Valiente-Banuet | Shade as a Cause of the Association Between the Cactus Neobuxbaumia Tetetzo and the Nurse Plant Mimosa Luisana in the Tehuacan Valley, Mexico[END_REF]: by using only point co-occurrence, this type of interaction is missed. On the other hand, other positive interactions such as protection against herbivory should lead to an increase in point co-occurrence only, because facilitated plants have to occur below the canopy of the facilitator.

Figure 5 . 1 :

 51 Figure5.1: Co-occurrence occurs where two plants are located at the same exact vertical point. A positive interaction can also give rise to two plants being close to each other more often than expected by chance -aggregation, yet never co-occur at the same position in space. In this case, however, the positive interaction cannot be inferred using only point co-occurrence.

5. 2 . 2 Figure 5 . 2 :

 2252 Figure5.2: The observed pattern of occurrence of a given species, at a given distance from another individual depends on multiple, superimposed mechanisms of interaction. Some of these mechanisms such as preemption of space only occur below the canopy of the individual (b), while others may have an effect that extends further (c, d). The net effect of these interactions is reflected in the probability of occurrence of the target species (e) -which can be measured on data.

  target in general. A more interaction-focused conservation could help us not miss the protection of this important, if not essential, facet of ecological communities. Part III P L AT E S & A P P E N D I C E S P L AT E S 123 Plate 1: General view of the landscape in La Crau, with a sheepfold in the background. Plate 2: Grazed vegetation in La Crau, with a sheep herd in the background. Plate 3: A large sheep flock gathering near the sheepfold in the evening (Peau de Meau). Plate 4: Typical Coussoul community, which is present in less-graed areas and dominated by the grass Brachypodium retusum. plates Plate 5: Close-up view of an individual of Onopordum illyricum displaying a protective effect on a diverse set of annual species growing under its canopy. Plate 6: Typical state of vegetation during the grazing season in front of a sheepfold (Peau de Meau): only few stress-tolerant species persist through the intense trampling. Plate 7: Ruderal species close to the sheepfold (La Grosse du Couchant) display a strong regrowth soon after the grazing pressure is released around sheepfolds (compare with Plate 6. Plate 8: The remaining foundations of a roman sheepfold (Vergières). The white line materializes the walls of the past sheepfold. plates Plate 9: Ruling tape used to measure the length over which plant individuals (or parts of them) intersect the transect. In this picture, the grass Brachypodium retusum and Carlina hispanica intersect the tape, interspersed with leaves of Dactylis glomerata. Here, a leaf of C. hispanica would be recorded as overlapping the transect from 44.0 cm to 46.5 cm (white arrows). Schematic representation of a pair of transects sampled in the field (Peau de Meau, Coussoul-type vegetation). Each rectangle materializes one individual with its height class, and the color represents the species (pebbles are in black). Full names of species are given in Appendix 1. A P P E N D I X 1 : S P E C I E S C O D E S U S E D I N T H I S W O R K

  ing as arid (average annual rainfall of 540 mm/year) as true steppe biomes. The maintenance of such a landscape is due to the combination if several local factors that reduce water availability and encroachment by ligneous species.

Figure 5 . 3 :

 53 Figure 5.3: General situation of the Crau. The main area is located within the dark triangle.

Figure 5

 5 Figure 5.4: (a) Aerial view of a sheepfold (Nouveau-Carton); (b) shematical illustration of the sheepfold within a grazing zone (b, continuous lines indicate roads). Sheep station at night in enclosures laid by the sheepfold (dark green patch in the aerial view), and driven within the zone during the day by the shepherd (b, dotted lines).

Figure 5 . 5 :

 55 Figure 5.5: Grazing pressure as a function of distance to the sheepfold in the S-E direction. Values are extracted from Dureau (1997) for his study of the "Coussouls de l'Opéra". With the exception of this study, the spatial distribution of grazing pressure within grazing areas remains (to my knowledge) undocumented. The fit is a function of the form y = a/(1 + b * x).

Figure 5 . 6 :

 56 Figure 5.6: Geographical distribution of the sampled sites in the Crau area.

Figure 5 . 7 :

 57 Figure 5.7: Distances of sampling for quadrats, transects and soil samples for each site.

Figure 5

 5 Figure 5.8: (a) Information recorded around each quadrat (p/a stands for presence/absence) (b) Spatial distribution of the 16 sub-samples for soil sampling.

Figure 5 . 9 :

 59 Figure5.9: Graphical summary of surveys carried out for a given site. Four vegetation belts have been identified based on a set of indicator species, as indicated by colors and labels. Labels indicate notable species present in each belt. Quadrats were carried out at different distances from the sheepfold entrance. In each belt, a pair of 5m transects was laid out to measure the spatial distribution of species individuals. Soil was also sampled at least four times at the average distance of each pair of transects.

  the center of the target individual of the interaction is under the canopy of the interactor) orβ i→j (d) = a ij * exp( r i -d δ ) if d > r i i.e. the center of the target individual is outside the canopy.

ParameterFigure 5 . 11 :

 511 Figure5.10: Geometric representation of interactions between plants. An interaction between two individuals is defined by a function β i→j between the pair of species that depend on a pairwise coefficient a ij , the distance between the two individuals d kl and δ that characterizes the distance of interaction. The interaction can increase locally the focal process (mortality or recruitment, red) or a decrease it (blue).

Figure 5 . 14 :

 514 Figure 5.13: Temporal variations of soil spatial change around sites

∆

  Figure 5.15: Graphical representation of the parameters used for priors.Here, y represents the variable under study (e. g. species richness) and d represents the distance to the sheepfold.

  capacité d'une communauté à répondre aux perturbations. Une première approche, basée sur l'usage d'indicateurs, a été utilisée pour identifier les communautés susceptibles de répondre brutalement à un changement de condition environmentale. Cependant, cette approche n'est applicable que lorsqu'un patron particulier est présent dans le système écologique -une boucle de rétroaction. Dans la plupart des cas, les interactions entre les espèces sont peu connues et il est donc nécessaire de les documenter afin de comprendre comment elles peuvent affecter la stabilité de la communauté. Depuis le début des années 2000 et le développement d'ouils permettant l'analyse de réseaux (Morales-Castilla et al., 2015; van Borkulo et al., 2015; Poisot et al., 2012; Delmas et al., 2017), de nombreuses études cherchent à cartographier les interactions entre espèces dans les communautés et de comprendre comment leur structure détermine la résilience du système en entier. Cette vision unidirectionnelle d'interactions données et fixes, déterminant la stabilité et la résilience d'une communauté écologique, a été cependant critiquée (de Ruiter, 2005). De fait, les interactions entre espèces sont variables, et affectées par les perturbations elles-mêmes.

Figure 5 . 16 :

 516 Figure 5.16: Un communauté écologique est définie par un ensemble d'espèces et leurs interactions. Ces interactions peuvent présenter certains patrons structurels (par ex. des boucles de rétroaction) qui détermine la résilience et la stabilité de la commmaunté -autrement dit sa réponse au perturbations. Cette dernière, combinée avec les caractéristiques effectives des d'une perturbation affecte en retour les abondances et les intensités d'interaction entre espèces

Figure 5 . 17 :

 517 Figure 5.17: La composition spécifique des communautés des zones sèches des clairières subalpine est bimodale (a), avec une communauté dominée par C. filifolia (2) contrastant avec une communauté plus diverse (1). Ces deux groupes diffèrent principalement en terme de granulométrie du sol (Gravierpebble -et sablesand) et de pourcentage de sol nu (Bare ground). La présence d'une bimodalité dans la composition spécifique à des niveaux hydriques similaires, et son lien avec la qualité du sol pourraient suggérer une boucle de rétroaction positive entre végétation et qualité du sol.

  étudiées présentent un déficit d'associations positives, et par conséquent un ratio entre associations positives et négatives toujours négatif. Cela suggère que la compétition est le mécanisme le plus structurant dans ces communautés de plante. Le stress biotique, et plus spécifiquement le pâturage et le piétinement par les moutons diminue la force et le nombre des associations entre espèces (Figure 5.18). Sous des niveaux importants de stress, les interactions entre espèces jouent un rôle vraisemblablement mineur sur la dynamique de la communauté végétale.

Figure 5 . 18 :Figure 5 Figure 5 . 20 :

 5185520 Figure5.18: Variations du ratio d'association (taille d'effet comparé à l'attendu aléatoire) en fonction de la distance à la bergerie. Le ratio est proche de l'attendu aléatoire proche de la bergerie, ce qui marque une diminution de l'importance des associations dans la structuration de la communauté dans ces zones à forte perturbation (pâturage très important). Ce ratio devient plus négatif à mesure que le pâturage diminue (la distance augmente), mettant en évidence une importance accrue de processus de compétition entre espèces.

Table 0 . 1

 01 Diversity in the measures of stability and resilience . . . . . . . . . . . . . . . . . . . . . . . Table 0.2 Structural patterns in interaction networs and community stability . . . . . . . . . . . . . . . . Table of the surveyed sites . . . . . . . . . . . . Parameters used in the model simulations . . .

	Table 5.2

5

Table

5

.1

Table 0 . 1 :

 01 Diversity in the measures of stability and resilience: a non-exhaustive list of stability metrics, with their usual domain of application in ecological work (data from[START_REF] Donohue | On the dimensionality of ecological stability[END_REF].

	Stability measure	Aspect quantified	Domain of most fre-	Note
			quent application
	Constancy (inverse of	Magnitude of variation over	Empirical work
	Variability)	time/space
	Persistence,	Appearance/extinction of	Variable	Robustness often quantifies the number of ex-
	Robustness	species over time.		tinctions after a species is removed from a sys-
				tem, i. e. the number of secondary extinctions
	Resistance	Capacity to remain	Empirical work
		unchanged under
		disturbance
	Asymptotic stability	Rate of return to equilibrium	Theoretical work
		after infinitesimal
		perturbation
	Engineering Resilience Speed of return to	Variable	This metric itself can be measured in different
		equilibrium after a realized		ways on empirical data, Ingrisch & Bahn (2018)
		perturbation
	Ecological Resilience	Maximum perturbation	Variable
		before a system exhibits
		qualitatively different
		functioning

  Continuous lines indicate stable states, and the dashed line indicates the limit of the two basin of attractions where multiple stable states are present for the same grazing values. This line delineates two basins of attraction, which are the set of values from which a system will converge to a particular stable state (grey arrows). White dots indicate critical points, at which one of the multiple stable states becomes unstable.

Table 0 .

 0 2: Some structural patterns found in interaction networks and their link with community stability or dynamics.

  [START_REF] Butitta | Spatial early warning signals in a lake manipulation[END_REF][START_REF] Berdugo | Plant spatial patterns identify alternative ecosystem multifunctionality states in global drylands[END_REF][START_REF] Cline | Early warnings of regime shifts: Evaluation of spatial indicators from a whole-ecosystem experiment[END_REF][START_REF] Ratajczak | Changes in spatial variance during a grassland to shrubland state transition[END_REF]. In the context of the currently available and upcoming spatially extensive data, spatial EWS have been suggested

[START_REF] Majumder | Inferring critical points of ecosystem transitions from spatial data[END_REF][START_REF] Souza | Forecasting Bifurcations from Large Perturbation Recoveries in Feedback Ecosystems[END_REF]

. Evaluation of EWS on empirical data has also proved to be challenging, and a number of studies have reported various limitations of EWS when evaluated on out-of-lab data, sometimes due to reasons such as a lack of sufficiently long and finely resolved data sets, or underlying stochasticity

[START_REF] Gsell | Evaluating early-warning indicators of critical transitions in natural aquatic ecosystems[END_REF][START_REF] Burthe | Do early warning indicators consistently predict nonlinear change in long-term ecological data[END_REF][START_REF] Chen | Rising Variability, Not Slowing Down, as a Leading Indicator of a Stochastically Driven Abrupt Transition in a Dryland Ecosystem[END_REF][START_REF] Weerman | Changes in diatom patch-size distribution and degradation in a spatially self-organized intertidal mudflat ecosystem[END_REF][START_REF] Moreno-De Las Heras | Assessing landscape structure and pattern fragmentation in semiarid ecosystems using patch-size distributions[END_REF][START_REF] Ashwin | Tipping points in open systems: Bifurcation, noise-induced and rate-dependent examples in the climate system[END_REF]

. Successful tests on ecosystem data remain rare -especially so for spatial EWS -which questions the usefulness of these indicators as tools to assess ecosystem degradation. The recent increase in the availability and extensiveness of remote sensing data provides a timely opportunity to test spatial EWS outside of controlled experiments (e.g.

[START_REF] Gsell | Evaluating early-warning indicators of critical transitions in natural aquatic ecosystems[END_REF][START_REF] Burthe | Do early warning indicators consistently predict nonlinear change in long-term ecological data[END_REF][START_REF] Eby | Alternative stable states and spatial indicators of critical slowing down along a spatial gradient in a savanna ecosystem[END_REF]

  ). Soil analyses included pH (in water), cation total exchange capacity (Metson method), total CaCO 3 , total organic matter, total phosphorus (Olsen method), exchangeable K, Mg, Ca and

	Belt A Ruderal species
		Malva sylvetris
		Chenopodium spp.
	5	Urtica urens
	Na and total Nitrogen.
		Belt B Thorny species
		Onopordum illyricum
		Carduus sp.
		Belt C Trifolium community
		Trifolium subterraneum
		Trifolium spp.
		Medicago spp.
	Grazing intensity	Belt D "Coussouls" community Brachypodium retusum
	increases	Asphodelus ayardii Sideritis romana
	N
	1x1m quadrat
	5m transects 2 per vegetation belt
	soil sampling (16 subsamples within a 5m radius)

  The β m 0 the β m s k →s j (d k ) represent the effect of individual k (of species s k ) on the mortality of individual j at distance

	d jk .
	recruitment of new individuals	At each time step, each
	individual spreads f propagules randomly in the whole system. Be-
	cause plants take physical space the recruitment process is divided in
	two steps: (i) physical access of the propagule to the ground and (ii)
	successful establishment of the new individual.

route des Junchas -26460 Bourdeaux, France

has been published in 2018 in Methods in Ecology and Evolution 1 , is under continuous development at its homepage 2 . 'spatialwarnings' is also widely distributed to R users through CRAN 3 .

The tools provided by this package were also used in modelling work, to investigate the role of different types of stressors on the predictive performance indicators. This has resulted in another publication in Ecological indicators 4 . https://doi.org/10.1111/2041-210X.13058 https://github.com/spatial-ews/spatialwarnings https://cran.r-project.org/package=spatialwarnings https://doi.org/10.1016/j.ecolind.2017.10.071

. found to detect more species and thus document more exhaustively the species composition of vegetation communities.

| |