Camille Saade
email: camille.saade@umontpellier.fr

Emanuel A Fronhofer

Benoit Pichon

Sonia Kéfi

Benoît Pichon

Landscape structure and dispersal rate drive large scale catastrophic shifts in spatially explicit metapopulations

Keywords: metapopulation, space, bistability, catastrophic shifts, hysteresis Figures: 6, Supplementary Figures: 12, Supplementary Tables: 1

Landscape structure and dispersal rate drive large scale catastrophic shifts in spatially explicit metapopulations

Introduction

As humans exert increasing pressures on ecosystems, understanding how they respond in terms of biomass, species diversity and composition is one of the most pressing issue in ecology. Importantly, some ecosystems can exhibit abrupt responses to gradual changes in environmental conditions [START_REF] Noy-Meir | Stability of grazing systems: an application of predator-prey graphs[END_REF][START_REF] May | Thresholds and breakpoints in ecosystems with a multiplicity of stable states[END_REF][START_REF] Scheffer | Catastrophic shifts in ecosystems[END_REF]. Such responses, referred to as "catastrophic shifts", are usually attributed to the existence of multiple stable ecosystem states between which ecosystems can shift following a perturbation or when a threshold in environmental condition, or "tipping point" is passed [START_REF] Suding | Alternative states and positive feedbacks in restoration ecology[END_REF]Kéfi et al., 2007b[START_REF] Kéfi | When can positive interactions cause alternative stable states in ecosystems?[END_REF]. Reverting the ecosystem to its initial state after a shift can be challenging if not impossible, a phenomenon known as "hysteresis". Catastrophic shifts have been described in various ecosystems all over the world [START_REF] Biggs | The regime shifts database: a framework for analyzing regime shifts in social-ecological systems[END_REF]. Common examples include the eutrophication of shallow lakes [START_REF] Scheffer | Alternative equilibria in shallow lakes[END_REF][START_REF] Meijer | Long-term responses to fish-stock reduction in small shallow lakes: interpretation of five-year results of four biomanipulation cases in the netherlands and denmark[END_REF][START_REF] Scheffer | On the dominance of filamentous cyanobacteria in shallow, turbid lakes[END_REF][START_REF] Carpenter | Management of eutrophication for lakes subject to potentially irreversible change[END_REF][START_REF] Jeppesen | Lake and catchment management in denmark[END_REF], the degradation of coral reefs [START_REF] Done | Phase shifts in coral reef communities and their ecological significance[END_REF][START_REF] Knowlton | Thresholds and multiple stable states in coral reef community dynamics[END_REF][START_REF] Mccook | Macroalgae, nutrients and phase shifts on coral reefs: scientific issues and management consequences for the great barrier reef[END_REF][START_REF] Nyström | Coral reef disturbance and resilience in a humandominated environment[END_REF], the transitions between woodlands and savannas [START_REF] Dublin | Elephants and fire as causes of multiple stable states in the serengeti-mara woodlands[END_REF][START_REF] Wilson | Positive-feedback switches in plant communities[END_REF][START_REF] Walker | Rangeland ecology: managing change in biodiversity[END_REF] and the desertification of grasslands [START_REF] Kassas | Desertification: a general review[END_REF]Rietkerk et al., 1997;[START_REF] Nicholson | Land surface processes and sahel climate[END_REF][START_REF] Wang | Ecosystem dynamics and the sahel drought[END_REF][START_REF] Reynolds | Global desertification: building a science for dryland development[END_REF]Kéfi et al., 2007a). Because they are difficult to predict and reverse [START_REF] Scheffer | Catastrophic regime shifts in ecosystems: linking theory to observation[END_REF] and because they can negatively affect human livelihoods [START_REF] Reynolds | Global desertification: building a science for dryland development[END_REF][START_REF] Biggs | The regime shifts database: a framework for analyzing regime shifts in social-ecological systems[END_REF], catastrophic shifts have gained a large interest in the literature, providing us with a theoretical framework rooted in bifurcation theory to describe the most common catastrophic shifts (e.g., desertification: [START_REF] Noy-Meir | Stability of grazing systems: an application of predator-prey graphs[END_REF][START_REF] May | Thresholds and breakpoints in ecosystems with a multiplicity of stable states[END_REF]Rietkerk and van de Koppel 1997;[START_REF] Klausmeier | Regular and irregular patterns in semiarid vegetation[END_REF]Kéfi et al. 2007b,a; or shallow lakes eutrophication [START_REF] Scheffer | Catastrophic shifts in ecosystems[END_REF][START_REF] Carpenter | Management of eutrophication for lakes subject to potentially irreversible change[END_REF]).

However, early theoretical work on catastrophic shifts mainly focused on isolated systems (but see [START_REF] Keitt | Allee effects, invasion pinning, and species' borders[END_REF][START_REF] Van Nes | Implications of spatial heterogeneity for catastrophic regime shifts in ecosystems[END_REF]. Real ecosystems are networks of connected entities between which matter and energy can flow [START_REF] Leibold | The metacommunity concept: a framework for multi-scale community ecology[END_REF]. This means that if one of these entities experiences a catastrophic shift locally, this shift can spread through the spatial network and possibly trigger other shifts. Scaling up our understanding of catastrophic shifts from the local 3 scale to the regional scale remains a challenge. In particular, whether local multistability implies regional multistability or whether different stable states can coexist in space (spatial multistability) and how this is modulated by a landscape's spatial structure are still open questions. Recent studies have started to address these questions by including space explicitly in models, under two main settings: i) in continuous space, which is well suited to relatively homogeneous habitats where clear spatial entities are not easily identified (see the example of Lake Veluwe in van de [START_REF] Van De Leemput | Resilience of alternative states in spatially extended ecosystems[END_REF], and ii) in discrete space, better representing habitats with clear spatial entities (or patches) connected through dispersal of individuals and flows of resources (e.g., an archipelago such as the Åland Islands in [START_REF] Hanski | Metapopulation persistence of an endangered butterfly in a fragmented landscape[END_REF]. As a consequence of these recent developments, we have a good understanding of how locally bistable dynamics can affect the regional behaviour of spatially structured ecosystems: alternative stable states are not expected to coexist in space [START_REF] Keitt | Allee effects, invasion pinning, and species' borders[END_REF][START_REF] Van De Leemput | Resilience of alternative states in spatially extended ecosystems[END_REF], unless dispersal is low and space is discrete [START_REF] Keitt | Allee effects, invasion pinning, and species' borders[END_REF], or strong stochasticity and heterogeneity are at play [START_REF] Van Nes | Implications of spatial heterogeneity for catastrophic regime shifts in ecosystems[END_REF][START_REF] Martín | Eluding catastrophic shifts[END_REF]. A spatially structured system will experience sharp transitions between fully occupied and fully empty states. These can usually not be called catastrophic shifts as regional bistability and hysteresis are largely reduced due to the dominance of a single of the two possible stable states [START_REF] Keitt | Allee effects, invasion pinning, and species' borders[END_REF][START_REF] Hilt | Abrupt regime shifts in space and time along rivers and connected lake systems[END_REF][START_REF] Van De Leemput | Resilience of alternative states in spatially extended ecosystems[END_REF]. However, all studies so far assume overly simplified spatial structures, either 1-D [START_REF] Keitt | Allee effects, invasion pinning, and species' borders[END_REF][START_REF] Van De Leemput | Resilience of alternative states in spatially extended ecosystems[END_REF] or 2-D continuous space [START_REF] Martín | Eluding catastrophic shifts[END_REF], linear [START_REF] Keitt | Allee effects, invasion pinning, and species' borders[END_REF][START_REF] Hilt | Abrupt regime shifts in space and time along rivers and connected lake systems[END_REF] or grid-like (van Nes and [START_REF] Van Nes | Implications of spatial heterogeneity for catastrophic regime shifts in ecosystems[END_REF] discrete systems. This implies that we currently don't know how more realistic spatial settings affect landscape-scale stability. Indeed, in real ecosystems, suitable habitat is usually neither continuous nor organized on a regular lattice but has a particular structure: terrestrial populations, for example, usually show emergent modularity [START_REF] Gilarranz | Generic emergence of modularity in spatial networks[END_REF] while riverine systems are typically dendritic [START_REF] Carraro | Generation and application of river network analogues for use in ecology and evolution[END_REF][START_REF] Rinaldo | River networks as ecological corridors: Species, populations, pathogens[END_REF]. Importantly, these specific spatial configurations have been shown to affect ecological dynamics [START_REF] Gilarranz | Effects of network modularity on the spread of perturbation impact in experimental metapopulations[END_REF][START_REF] Altermatt | Dispersal in dendritic networks: Ecological consequences on the spatial distribution of population densities[END_REF][START_REF] Carrara | Dendritic connectivity controls biodiversity patterns in experimental metacommunities[END_REF]. This omission of spatial complexity is an important shortcoming of the current state of the literature as the properties of such habitats may 4 change how local bistability affects regional scale dynamics and equilibria [START_REF] Van Nes | Implications of spatial heterogeneity for catastrophic regime shifts in ecosystems[END_REF][START_REF] Van De Leemput | Resilience of alternative states in spatially extended ecosystems[END_REF].

Here, we aim to fill this gap by studying how local bistability affects the regional scale dynamics of spatially complex populations, using a metapopulation framework in which local patches are connected to each other via dispersal of individuals. We simulate the dynamics of metapopulations with various spatial structures, ranging from classical abstract structures (lattice and random, Erdős-Rényi, spatial networks) to structures rooted in real systems (Random Geometric Graphs for terrestrial systems and Optimal Channel Networks for riverine systems; [START_REF] Gilarranz | Generic emergence of modularity in spatial networks[END_REF][START_REF] Carraro | Generation and application of river network analogues for use in ecology and evolution[END_REF]. For a range of dispersal rates, we measure the size of hysteresis and the position of tipping points at the scale of the whole landscape. We also study how the structure of perturbations themselves affects landscape dynamics by conducting targeted local perturbations of firstly highdegree patches vs. low degree patches and secondly neighboring patches vs. independent patches.

Model description

Metapopulation model

We used a model of a spatially structured metapopulation describing the dynamics of n patches linked by the dispersal of individuals. Here is a short description of the model, but see the supplement S1 for a more in-depth presentation. The dynamics of the vector -→ x (containing the local biomasses x i ) are described by a system of n coupled differential equations, which can be written in a matricial form:

d - → x dt = r - → x 1 - - → x K -B - → x A + - → x + µ M - → x (1)
The first two terms of Eq. 1 describe the local dynamics of the biomass, while the last term describes the dispersal between patches.

Local dynamics

We describe the local dynamics using a model derived by [START_REF] Noy-Meir | Stability of grazing systems: an application of predator-prey graphs[END_REF] in which plants grow logistically with a growth rate r and carrying capacity K, and are harvested 5 by a fixed consumer with a Holling type II functional response (B: maximal harvesting rate, A: half saturation constant). For the rest of this work, we used the same fixed values for r, K and A (Table S6.1). These parameter values give rise to two tipping points at B = 50 and B = 56.25 between which the system is bistable: it admits two stable equilibria (x i = 0 and

x i = r 1 > 0)
separated by an unstable equilibrium (0 < r 2 < r 1 ) (Fig. 1a).

Note that this is one of several possible models of bistability that we chose because of its simplicity and intuitive biological interpretation. Importantly, we show in the supplement S3 how this model and three other common models of bistable systems share a common structure through parameter aggregation which highlights the generality of our findings.

Spatial dynamics

We consider that individuals leave their local patch at a constant rate µ and that the dispersing individuals are distributed evenly among adjacent patches. Mathematically, the vector describing net dispersal for all patches is given by the last term of Eq. 1 (µM -→ x ) where M is a matrix describing the structure of the landscape: the diagonal terms describe emigration (M i,i = -1)

and the off-diagonal terms describe the neighbourhood relationship between patches (see detail in supplement S1).

Integration scheme We simulated the dynamics using the function ode (R-package deSolve, version 1.28). We performed linear regressions on the biomass of each patch over the last 100 simulated values to determine whether the metapopulation had reached an equilibrium (slope smaller than 10 -3 ).

Analysis of a two-patch system

Before moving on to large systems which are only accessible through simulation (see section "Landscape design"), we conducted analytically the stability analysis of a simple system for demonstrating general principles. We used a particular case of Eq. 1 with only two patches exchanging biomass at a rate µ (supplement S1: "Two-patch system"). We determined all its existing equilibria as well as 6 their stability using MatCont (version 7.2, in MATLAB R2017b).

Landscape design

We used four different types of landscapes (see supplement S2 for more details on landscape generation, and Fig. S6.1 to S6.3 for examples of landscapes): i) Regular graphs, where patches are arranged regularly on a lattice with periodic boundary conditions. While not very realistic, such landscapes have been used in previous studies [START_REF] Keitt | Allee effects, invasion pinning, and species' borders[END_REF][START_REF] Van Nes | Implications of spatial heterogeneity for catastrophic regime shifts in ecosystems[END_REF] and we include them for the sake of comparison. ii) Erdös-Rényi graphs are obtained by randomly connecting pairs of patches with a fixed probability p and can be thought of as a null model. iii) Random Geometric Graphs (RGGs) are obtained by randomly drawing patch coordinates in space and connecting patches depending on distance. The random distribution of patches in space results in the emergence of modularity which is thought to be characteristic of terrestrial systems [START_REF] Gilarranz | Generic emergence of modularity in spatial networks[END_REF]. iv) Optimal Channel Networks (OCNs) are generated by simulating geomorphological processes to obtain a structure closely resembling that of a river [START_REF] Carraro | Generation and application of river network analogues for use in ecology and evolution[END_REF].

For each type of landscape, we generated 50 replicates with each n = 100 patches. The main text analysis was made on landscapes where patches had on average approximately 4 neighbours (except for OCNs as their generating process constrains their connectivity), and we conducted sensitivity analysis on the connectivity using networks with higher and lower connectivity (supplement S4).

Simulations

For each type of network, we generated 50 replicates and computed their bifurcation diagram for 100 values of dispersal rates (µ) between 0.001 and 1. We used 46 values of harvesting rates (B) equally spaced between 50 and 56.25 (the bistability range of a single patch). For each value of B, we computed the high-biomass (resp. low-biomass) branch of the bifurcation diagram by initially setting each patch to the positive (resp. null) equilibrium of a single patch. Once the system had reached an equilibrium, we degraded (resp. restored) 5% of the patches by setting them to a null 7 biomass (resp. to their positive equilibrium). We then measured the steady-state average biomass in order to plot the bifurcation diagram.

We used three different perturbation schemes: Firstly, we chose the 5 patches to be degraded/restored randomly (Fig. 234). Secondly, we focused on the influence of node degree by degrading/restoring patches chosen from either the most or least connected nodes (Fig. 5). We conducted these perturbations in Erdös-Rényi graphs, RGGs and OCNs but not in regular graphs (as their nodes all have the same degree). Thirdly, we focused on the influence of the spatial structure of perturbations by degrading/restoring patches that were either clustered or dispersed in space (Fig. 6). In regular graphs, this was done by drawing the 5 perturbed nodes either from a 9 patch neighborhood (3x3 patches in grid landscapes or 9 adjacent patches in circular landscapes) for the "clustered" modality or drawing them from the whole landscape for the "dispersed" modality. In

RGGs and OCN, we identified modules using the function edge.betweenness.community (R-package igraph version 1.2.6) and the dispersal matrix M, and drew 5 patches from a single modules for the "clustered" modality or from 5 different modules for the "dispersed" modality. We did not conduct this third analysis on Erdös-Rényi graphs as they are not spatially explicit so there is no notion of proximity between patches.

Characterizing large scale shifts

From the established bifurcation diagrams, we extracted synthetic information on the large scale behaviour of the system.

State diagram

We affected one of three categorical states to each simulation depending on its steady-state average biomass: it was considered "fully occupied" when the average biomass was higher than 99% of the pre-perturbation high-biomass equilibrium, "empty" when the average biomass was less than 1% of it and "partially occupied" in between. We constructed state diagrams for each type of networks by averaging the states borders over all replicates and plotting them as a function of the harvesting rate (B) and dispersal rate (µ).

8

Hysteresis size We computed the size of the hysteresis as the area in between the higher and lower branches of the bifurcations diagram.

Tipping point Finally, we computed the position of the tipping points from the bifurcation diagrams as the value of harvesting rate (B) at which the average biomass decreased (for the degradation trajectory) or increased (for the restoration trajectory) the most.

Results

Analysis of a two-patch system

Starting from an isolated patch exhibiting bistability (Fig. 1a) and connecting it to another patch, we found that, at high dispersal, spatial bistability -the stable coexistence of patches in different states -is not possible and the whole system behaves as a single bistable unit (Fig. 1b, supplement S5).

As dispersal decreases, spatial bistability becomes possible as new stable states appear (Fig. 1c, supplement S5) with strong differences in biomass between the two patches.

Stable states of larger metapopulations

We now focus on larger metapopulations made of 100 patches. The results described in the following sections held qualitatively when considering graphs that were overall more or less connected (see sensitivity analysis in supplement S4).

On a regular grid, a larger metapopulation can reach different states (fully occupied, empty or partially occupied) depending on dispersal (µ) and harvesting (B) rates (Fig. 2). At high dispersal rates (Fig. 2a-d, µ > 0.15), only two stable states are possible: the landscape is either fully occupied (the average biomass is that of a single patch) or empty (the average biomass is null), and the metapopulation is bistable as a whole (regional bistability). The bifurcation diagram of the average biomass is then qualitatively similar to that of a single patch (Fig. 2d), albeit with a displaced restoration tipping point. As dispersal rate decreases, the bistability domain shrinks and partially occupied states (where some of the landscape patches are empty and others are occupied) become stable (i.e., spatial bistability; Fig. 2c,e). At very low dispersal rates, the system behaves roughly as a collection of independent patches: spatial bistability is common, and local degradation (resp.

restoration) can neither spread in space nor be reversed. As a consequence, the bifurcation diagram of the average biomass mostly reflects the initial conditions (Fig. 2f).

These results hold qualitatively for other types of landscapes such as Erdös-Rényi and Random Geometric Graphs (RGG), even though spatial bistability is then more common than for regular grids (Fig. 3a-c). In Optimal Channel Networks (OCNs), we observed a similar state space diagram, although the bistability domain was reduced and did not expand much with increasing dispersal rates, and the partially occupied (i.e., spatially bistable) state was stable at higher dispersal rates than for the other networks (up to µ ≈ 0.25) (Fig. 3d).

Characterisation of regional scale catastrophic shifts

Hysteresis size For all landscape types, the hysteresis size of the metapopulation was smaller than that of an isolated patch (Fig. 4a). Hysteresis sizes were similar across network types at low dispersal rates (µ < 0.1) but quickly diverged for higher dispersal rates: regular graphs had the larger The colors denote the final state of the system: "fully occupied" (green) means that the total equilibrium biomass was greater than 99% of the maximal equilibrium biomass, "empty" means that the total equilibrium biomass was smaller than 1% of the maximal equilibrium biomass, and "partially occupied" means that the system was neither empty nor fully occupied. optimal chanel networks (OCNs). Each panel was obtained by combining the state space diagram of a landscape starting with 5% of random occupied patches and the state space diagram of a landscape starting with 5% of random empty patches, averaged over 50 replicates, similarly as Fig. 2(a-c).

with the dispersal rate to reach a plateau (µ > 0.1). Hysteresis size in OCNs also decreased with dispersal rate at first, but stayed at a constant value for higher dispersal rates (µ > 0.2).

Tipping points For all landscape types, the degradation (resp. restoration) tipping points happened at a lower (resp. higher) harvesting rate than in an isolated patch, which explains the lower hysteresis size in spatially structured landscapes (Fig. 4b andc). The degradation tipping point was mostly unaffected by dispersal rate beyond a threshold (µ > 0.1 for regular graphs, Erdös-Rényi graphs and RGGs and µ > 0.25 for OCNs). The restoration tipping point, however, was more affected by dispersal rate (Fig. 4c). Interestingly, this means that the patterns observed for hysteresis size are mainly explained by modifications of the restoration rather than of the degradation point. OCNs were the landscapes whose tipping points were overall the most different from an isolated patch.

Similarly to what was observed for hysteresis size, the tipping points were closer to those of a single patch at low dispersal rates, but then deviated from it before reaching a plateau.

Perturbation modality and large scale catastrophic shifts

Lastly, we considered how different types of perturbations affected large scale shifts, as opposed to the random perturbations considered up to this point.

Spatial proximity of perturbations Both in regular graphs and RGGs, spatially clustered perturbations resulted in a smaller hysteresis compared to spatially dispersed perturbations (Fig. 5a) for all dispersal rates higher than ∼ 0.1. This was once again explained by differences in the restoration tipping point: while degradation happened at the same harvesting rate in both perturbation modalities (Fig. 5b), the restoration happened at higher harvesting rates for clustered perturbations 13 (Fig. 5c).

In OCNs, the spatial organization of perturbations had no effect on the hysteresis size (Fig. 5a) as both degradation and restoration were mainly unaffected by the perturbation modality (Fig. 5b andc).

Connectivity of perturbed patches

In RGG, perturbing highly vs. lowly connected patches had no discernible effect on the hysteresis size and tipping points of the metapopulation (Fig. 6). In Erdös-Rényi graphs and OCNs, perturbing lowly connected patches resulted in a slightly smaller hysteresis at high dispersal rates (µ > 0.1) (Fig. 6a), once again because of differences in the restoration point: restoration happened at higher harvesting rate for perturbations of lowly connected nodes compared to highly connected nodes (Fig. 6c), while the degradation point was unaffected (Fig. 6b). 14

Discussion

Stable states of metapopulations composed of locally bistable patches

We found that the stable states of metapopulations made of locally bistable patches depended strongly on dispersal rate and were qualitatively similar for various spatial structures and network connectivity.

Strong dispersal homogenized biomass over space and drove the system towards one of two stable states: either all the patches were occupied and the average biomass was roughly the same as the positive equilibrium of an isolated patch, or all patches were empty and the average biomass was null. These two states had overlapping stability ranges, making the metapopulation bistable as a whole, with a clear hysteresis and abrupt shifts. Therefore, we expect the local degradation or restoration of a few patches in an otherwise homogeneous metapopulation to be quickly reversed, which means that, in a metapopulation with strong dispersal, local restoration efforts are bound to fail unless the environmental conditions cross a metapopulation-scale tipping point. Once this tipping point is crossed, local restoration efforts should spread in space until the whole metapopulation is 15 restored.

When dispersal was sufficiently weak, patches were functionally independent from each other.

The exchange of biomass was too low for patches to impact the state of their neighbours and spatial bistability was common as occupied and empty patches could thus coexist in space. In our setting, perturbing only 5% of the patches, this resulted in a large scale behaviour very similar to that of a single patch: the unperturbed patches (95% of the landscape) shifted only when the tipping point of an isolated patch was reached, independently of the perturbations. From a conservation standpoint, the effective independence of patches ensures that a local degradation won't spread in space and policies could focus on local measures of protection without worrying about spatial effects. However, in the absence of significant fluxes between patches, spatial heterogeneity and stochasticity should have a strong impact on local dynamics and extinction probability. How this can affect large scale dynamics is still an open question as these processes can synergize in surprising ways: for example, [START_REF] Martín | Eluding catastrophic shifts[END_REF] found that demographic stochasticity coupled with a low dispersal rate results in smooth transition in 1-and 2-dimensional continuous space.

Interestingly, this means that in metapopulations made of locally bistable patches, strong dispersal should raise concerns as it translates into possible regional scale catastrophic shifts with a pronounced hysteresis and little prospect for local restoration. On the other hand, weak dispersal should open the door to local conservation and restoration efforts. This should be taken into account when planning conservation measures such as assisted migration or ecological corridors, as increasing a species dispersal rate could make it prone to large scale catastrophic shifts.

Characteristics of regional scale catastrophic shifts

We observed a smaller hysteresis size than expected from the dynamics of isolated patches for all landscape types. This is particularly pronounced in OCNs where hysteresis was largely reduced (less than 1/10 th of the hysteresis of a single patch for dispersal rates > 0.2). This was expected from previous studies, which predicted that hysteresis should mostly disappear at high dispersal in one dimensional continuous space (van de [START_REF] Van De Leemput | Resilience of alternative states in spatially extended ecosystems[END_REF] or in linear metapopulations (Keitt 16 et al., 2001;[START_REF] Hilt | Abrupt regime shifts in space and time along rivers and connected lake systems[END_REF]. Since OCNs are made of large linear parts, they fit this prediction best. The other types of landscapes have more links per patch on average than OCNs, which diluted the effect of any focal patch on more than two neighbours. As a consequence, a local degradation or restoration in an otherwise homogeneous landscape is less likely to spread to its neighbours and more likely to be reversed than in less connected landscapes, resulting in a larger hysteresis.

Interestingly, this reduced hysteresis compared to an isolated patch is mostly due to an earlier restoration, i.e., restoration was reached at harsher conditions than expected from the local dynamics.

This is due to an asymmetry in the bifurcation diagram of local dynamics: the tipping point towards restoration is a transcritical bifurcation (Fig. 1a, point F1), while the tipping point towards degradation is a saddle-node bifurcation (Fig. 1a, point F2). As a consequence, the null biomass equilibrium has a small basin of attraction while the positive equilibrium has a large basin of attraction, even close to its tipping point, so local restorations are more easily spread in space than local degradations. This asymmetry of the basins of attraction holds true for various models of desertification (Rietkerk and van de Koppel, 1997;[START_REF] Klausmeier | Regular and irregular patterns in semiarid vegetation[END_REF]; we therefore expect these systems to also exhibit a reduced hysteresis and easier restoration in metapopulations than expected from the dynamics of a single patch. However, it is important to keep in mind that other models of bistability should behave differently. For example, a system with two transcritical bifurcations should show almost no hysteresis, while a system with saddle-node bifurcations should show an hysteresis of similar size to an isolated patch. Interestingly, a system with a transcritical bifurcation for the degradation and a saddle-node bifurcation for the restoration should have an hysteresis size similar to what we find here, but mainly through an earlier degradation.

Hysteresis size was found to vary with dispersal rate. The largest hysteresis size was observed at the lowest dispersal rate (µ = 0.001) for all types of networks, and quickly decreased with dispersal rates. After this initial decrease, the hysteresis size increased with dispersal to a plateau. This means that the smallest hysteresis size and the easiest restoration happened for intermediate dispersal rates Characterizing how hysteresis size varies with dispersal rate is important from a conservation standpoint, because a smaller hysteresis means that it is easier to restore the system after a large scale degradation. We can try to derive broad predictions about which species should raise concern for large scale catastrophic shifts. In terrestrial systems, dispersal increases with height for plants [START_REF] Tamme | Predicting species' maximum dispersal distances from simple plant traits[END_REF] and often with body size for animals [START_REF] Stevens | A comparative analysis of dispersal syndromes in terrestrial and semi-terrestrial animals[END_REF], while population growth rate is expected to decrease with body mass [START_REF] Savage | Effects of body size and temperature on population growth[END_REF]. Because of these relationships, we can speculate that one should be more concerned by the possibility of large scale catastrophic shifts in systems dominated by large species. Similarly, oceanic systems may be more prone to large scale catastrophic shifts with a marked hysteresis since oceanic species can exhibit strong propagule dispersal that impacts local dynamics [START_REF] Kinlan | Propagule dispersal in marine and terrestrial environments: a community perspective[END_REF].

(µ ≈ 0.1)
Of course, dispersal also depends strongly on other species traits (e.g., seed dispersal syndromes, body temperature) and such broad predictions remain speculations and cannot replace the detailed knowledge on a species of interest.

The landscape type also affected the metapopulation regional scale behaviour. In OCNs, hysteresis size did not increase with dispersal rates but stayed low after the initial decrease. This is because OCNs are almost one-dimensional, with large linear parts. This linearity limited the aforementioned dilution of local restoration efforts and resulted in tipping points almost unchanged 18 by dispersal at high dispersal. On top of that, patches in OCNs have a low average degree, with most nodes having only 1 or 2 neighbours. Because of this, the dispersal from any focal patch was concentrated on a few neighbours. This allowed an easier induction of regional scale catastrophic shifts through a domino effect (even for degradation), as evidenced by the degradation and restoration tipping points happening at respectively lower and higher harvesting intensity than in the other landscape types. Interestingly, two other studies support this result: Hilt et al. ( 2011) modeled a river as a strictly linear chain of discrete patches and found that hysteresis disappears entirely at high dispersal rates, while [START_REF] Martín | Eluding catastrophic shifts[END_REF] showed that one dimensional systems tend to have smaller hysteresis and smoother transitions than their two (and even three) dimensional counterparts.

The other landscape types (grid and RGGs) have no linear parts and an higher average degree (4 neighbours on average). Because of this, the effect of local perturbations was split between several neighbours and diluted in the metapopulation, making it more difficult to induce regional scale shifts and resulting in a larger hysteresis. This was supported by our sensitivity analysis on network connectivity (supplement S4): networks with a smaller average degree (Fig. S4.1) had a smaller hysteresis while networks with a larger average degree (Fig. S4.2) had a larger hysteresis.

From these results, we expect large-scale catastrophic shifts with a pronounced hysteresis to be more likely in well-connected and two dimensional systems (e.g., terrestrial systems, open ocean) than in systems that are lowly connected or restricted to large linear stretches (e.g., rivers, coastal ecosystems).

Lastly, we explored how the spatial structure of perturbations themselves could affect large scale catastrophic shifts. We found that close-by restoration efforts, e.g., patches from a single module in RGGs, are more likely to induce regional scale restorations and results in a smaller hysteresis size than restoration efforts conducted randomly in space. This is because restoring close-by patches concentrates the added biomass on a smaller area (e.g., the module where we restore patches in RGGs) and limits the dilution of biomass in space. This makes it easier to trigger the large scale restoration of the metapopulation. Studies in continuous space highlight a similar result: van de Leemput et al. ( 2015) identifies a "minimal size of disturbance to initiate a travelling wave". Several disturbances smaller than this size won't initiate a large scale shift as they get diluted individually, but grouping these disturbances over a single area will initiate a large scale shift.

On the other hand, the connectivity of perturbed patches slightly increased the hysteresis size in Erdös-Rényi graphs and OCNs, and had surprisingly no effect in RGGs. This is probably because sampling from low degree patches for the perturbations incidentally results in sampling far away patches, which partially masks the effect of patch connectivity on hysteresis.

Conclusion

In conclusion, we showed that one can not extrapolate the bifurcation diagram of a single patch to predict regional scale dynamics. Instead, we highlight that a metapopulation can exhibit large scale catastrophic shifts and hysteresis, and that the position of its tipping points depends on both its spatial architecture and its dispersal rate. We find that restoration after a large scale degradation is easier with i) local restoration efforts concentrated on a restricted area and ii) intermediate dispersal rates. Our findings differ markedly from the predictions established on regular one-dimensional systems and show that we should consider the explicit structure of metapopulations or at least their properties (dimension, average degree and degree heterogeneity) when trying to predict their regional scale dynamics.

20

Spatial dynamics We consider that the individuals leave the patch i at a constant rate µ and are equally split between its neighbours. The patch also receives individuals from its neighbours in the same pattern which yields the following equation:

g(x i , x • ) = -µx i + µ ∑ j ad jacent i x j Ne j
where Ne j is the number of neighbours of j.

Global dynamics of the system The dynamics of the system are thus described by a system of n coupled differential equations, which can also be written in a matricial form:

d - → x dt = r - → x 1 - - → x K -B - → x A + - → x + µ M - → x (4)
where -→ x is the vector of local densities and M is a matrix describing the structure of the population:

the diagonal terms describe the emigration (M ii = -1) and the off-diagonal terms the immigration (M i j,i̸ = j = 0 if i and j are not neighbors and M i j,i̸ = j = 1/Ne j if i and j are neighbors).

Two-patch system

The two patch system we used for the analytical stability analysis of a simple spatial system is a particular case of the model presented above, described by the following equations:

dx 1 dt = rx 1 1 - x 1 K - Bx 1 A + x 1 -µ(x 1 -x 2 ) dx 2 dt = rx 2 1 - x 2 K - Bx 2 A + x 2 -µ(x 2 -x 1 )

S2 Landscape design

We generated landscapes of n = 100 patches. We used four different landscape structures in order to explore the effect of space on catastrophic shifts: regular graphs (e.g., patches arranged either linearly or on a lattice), Erdős-Rényi graphs (where patches are connected randomly), random geometric graphs (RGG, where patches are drawn randomly in space and connected depending on their distance) and optimal channel networks (that are generated using geomorphological processes).

The main text analysis is made on networks where the patches have on average 4 neighbors each (when possible), and we conducted sensitivity analysis on the connectivity by generating networks with higher and lower connectivity (section S4).

Regular graphs Regular graphs are graphs in which all nodes have the same number of neighbours: patches are either organized linearly on a cirlce (in 1-D) or on a regular lattice on a torus (for 2-D system). These graph have been used extensively for modelling explicit space [START_REF] Keitt | Allee effects, invasion pinning, and species' borders[END_REF][START_REF] Van Nes | Implications of spatial heterogeneity for catastrophic regime shifts in ecosystems[END_REF], so we chose to use them for comparison with existing studies.

They can be thought of as a null model of the most ordered spatial system possible. We used patches on a 10*10 lattice with periodic boundary conditions and connections to the 4 nearest neighbours for the main text. For the sensitivity analysis, we used patches on a 10*10 lattice with periodic boundary conditions and connections to the 8 nearest neighbours (2-D system, mean degree = 8)

and 100 patches on a circle with connections to the 2 nearest neighbours (1-D system, mean degree = 2).

Erdős-Rényi graphs Erdős-Rényi graphs are a class of graphs obtained by randomly connecting nodes: each pair of patches is connected with a probability p. The expected number of links per node (mean degree) is thus p * (n -1) where n is the total number of patches. We used Erdős-Rényi graphs with n = 100 patches and a connection probability of p = 4/(n -1) for the main text (expected mean degree = 4, realized mean degree = 4). For the sensitivity analysis, we used connection probabilities of p = 8/(n -1) (expected mean degree = 8, realized mean degree = 8) 28 and p = 2/(n -1) (expected mean degree = 2, realized mean degree = 2.7). After generating a graph, we checked if it was connected (i.e., that there existed at least an indirect path between all pairs of nodes) and discarded the graphs that were not connected.

Random geometric graphs Random geometric graphs (RGG) are obtained by drawing nodes randomly in space and connecting the pair of nodes that are closer than chosen a threshold distance.

We generated each RGG by drawing 100 pairs of x and y coordinates uniformly in [0, 1]. We then determined the threshold distance that yielded the average degree the closest to a target degree (main text: target degree = 4; sensitivity analysis: target degree = 8). We once again checked if the graphs were connected and discarded those that were not connected. Note that the connectivity of a RGG is constrained by the initial drawing of coordinates, so it was not always possible to reach the targeted average degree. RGGs connectivity deviated a bit from the regular graph counterpart (main text: realized average degree = 4.7; realized degree = 7.9)

Optimal channel networks Optimal channel networks are obtained by simulating geomorphological processes and in order to capture the structural properties of riverine networks. We generated OCN of n = 100 patches using the R-package OCNet (version 0.4.0) [START_REF] Carraro | Generation and application of river network analogues for use in ecology and evolution[END_REF]. Because the structure of OCN is constrained by the underlying generating process and generally scale-invariant, we could not manipulate the average degree of these networks (average degree = 1.98).

29

S3 Polynomial form of three classical desertification models

In this annex, we highlight the mathematical similarities between four classical models of bistable systems -three models of desertification taken from [START_REF] Noy-Meir | Stability of grazing systems: an application of predator-prey graphs[END_REF], [START_REF] Klausmeier | Regular and irregular patterns in semiarid vegetation[END_REF] and van de Koppel et al. (1997); and a model of logistic growth with an Allee effect. We argue that, despite the fact that they are built on different mechanistic processes, they behave similarly as they are all build from a simple third degree polynomial. In the first part, we conduct the exact stability analysis of this third degree polynomial. We then show in the subsequent parts how to relate models of desertification to this polynomial, and how it simplifies their stability analysis.

S3.1 Stability analysis of a dynamical model described by a third degree polynomial

The most straightforward way to model bistability is using a cubic polynomial to describe the evolution of a quantity (number of individuals, plant cover, share of users adopting a new technology...) in a closed system. Let us call P such a polynomial:

P : R → R x → dx dt = α 3 x 3 + α 2 x 2 + α 1 x + α 0
Realistically, α 3 should be negative (otherwise dx/dt → ∞ when x is high, meaning the population is unbounded). Hence, α 3 is restricted to negative values in what follows. In biological systems, realism also imposes a restriction on α 0 : since the system is closed and x describes a density of organisms that cannot be spontaneously generated, P(0) must be null, and hence α 0 = 0.

Additionally, x can only be positive when describing a population, so P only needs to be defined on

R + .
Hereafter, P refers to the previously defined polynomial with these additional restrictions:

P : R + → R x → dx dt = α 3 x 3 + α 2 x 2 + α 1 x, α 3 < 0
Existence of alternative equilibria: Since α 0 = 0, the system has an immediate equilibrium for x = 0. The other equilibrium can thus be found by solving:

α 3 x 2 + α 2 x + α 1 = 0 This is a degree two polynomial with discriminant ∆ = α 2 2 -4α 3 α 1 . If ∆ < 0 (i.e. 4α 3 α 1 > α 2 2 ),
the polynomial has no real root and 0 is the only equilibrium of this system (Fig. S3.1c)

If ∆ > 0, then there are two other roots to the system (that are only relevant when they are positive):

r 1 = -α 2 - √ ∆ 2α 3 r 2 = -α 2 + √ ∆ 2α 3
The nature of the system depends on the sign of these roots: if both are negative, then 0 is a stable equilibrium and the only biologically significant equilibrium. If both are positive, 0 is still a stable equilibrium but there is also a stable positive equilibrium (r 1 ) separated from 0 by an unstable equilibrium (r 2 ) and thus the system is bistable (Fig. S3.1d). If only one of them is positive, then 0 is unstable and r 1 is stable and the only positive equilibrium (Fig. S3.1e).

31

Sign of r 1 : Since α 3 is negative, r 1 is of the same sign as α 2 + √ ∆. Thus, r 1 is positive iff

α 2 > - √ ∆ ⇔                                α 2 > 0 or                α 2 < 0 and 4α 3 α 1 < 0
Thus r 1 is only negative when α 2 < 0 and 4α 3 α 1 > 0, i.e. in the top left corner of the α 2 -4α 3 α 1 plane (Fig. S3.1a).

Sign of r 2 : Similarly, the sign of r 2 is that of

α 2 - √ ∆. r 2 is positive iff α 2 > √ ∆ ⇔                α 2 > 0 and 4α 3 α 1 > 0
r 2 is thus only positive in the top right corner of the α 2 -4α 3 α 1 plane (FigS3.1a).

These results are summed up in Fig. S3.1:

• if α 1 is positive (crossed area in Fig. S3.1a), the system always reach a stable positive equilibrium (r 1 , Fig. S3.1g)

• If α 1 is negative, the behaviour of the system depends on the relative values of α 2 and 4α 1 α 3 : if α 2 2 > 4α 1 α 3 (white area in Fig. S3.1a), the system is bistable and can either reach 0 or a positive equilibrium (r 1 ) depending on its initial state (Fig. S3.1f). If α 2 2 < 4α 1 α 3 , the system has no strictly positive equilibrium and necessarily goes to 0.

32

Stability analysis and interpretation: One way to make light of this is to see that α 1 drives the system behaviour when x is close to 0 while α 2 and α 3 become more predominant at high x value, thus:

• if α 1 is positive, P is positive when x is close to 0 and thus 0 is unstable. Since P is negative when x is high (because α 3 is negative), P necessarily becomes null once (and only once) for a x > 0 which is a stable equilibrium.

• If α 1 is negative, P is negative around 0 and thus 0 is a stable equilibrium. In that case, either α 2 makes a positive contribution high enough (meaning α 2 2 > 4α 1 α 3 ) for P to reach 0 two more times (because P ultimately decreases to -infinity at high x) and the system is bistable or α 2 contribution is not enough (α 2 2 < 4α 1 α 3 ) to bring P to 0 again and thus 0 is the only equilibrium.

An intuitive way to sum these relationships is by looking at the α 2 -4α 1 α 3 plane (Fig. S3.1a): in the area under the x-axis (in green), the system necessarily goes to a positive equilibrium. Over the x-axis, the curve defined by y = x 2 separates two areas: on the left (in red) where 0 is the only stable equilibrium and the system goes to extinction. On the right (in orange), the system is bistable. By increasing α 2 (moving from left to right) or α 1 (moving from top to bottom), we drive the system away from extinction and towards either a positively stable domain (in green) or a bistable domain (in orange) respectively.

S3.2 Application to biological systems:

While some theoretical studies chose to use a third degree polynomial as a simple way to obtain bistable dynamics, models rooted in biological mechanisms are rarely expressed as third degree polynomials -because cubic terms cubic terms are hard to intuitively link to a process. Instead the dynamics of biological systems (dx/dt, usually expressing biomass or species density) are often expressed as a combination of: In the red area, 0 is the only equilibrium in R + and the system always goes to extinction (see panel c).

(c) (d) (e)
In the orange area, the system has three equilibrium in R + ([0, r 1 , r 2 ]) and goes either to extinction or r 2 depending on the initial x (see panel d). In the green area, 0 is an unstable equilibrium and the system always goes to r 1 that is stable (see panel e). The full line shows a transect happening when increasing both α 2 from 0.8 to 2.4 and α 1 from -0.6 to 0.2 (α 3 is fixed at -0. • linear terms (λ * x) to describe density independent processes (e.g., death rate in birth-death processes, growth in an exponential model).

• square terms to describe density dependent processes such as mass action law (e.g., intraspecific competition in a logistic growth model: r(xx 2 /K)).

• a quotient of low order polynomials (1 or 2) to describe saturating processes (e.g., type 2 functional response: n/(1 + n)).

We found that a number of these models can be restated as a quotient of two polynomials:

f (n) := dx dt = P(x) g(x)
where P is the third degree polynomial described earlier and g is a strictly positive and monotonously increasing second order (at most) polynomial (defined on R + ). Because of these properties, P and f share the same sign and roots despite f not being a polynomial. In the following section, we illustrate this using three classical models of desertification as well as a model of logistic growth with an Allee effect. We argue that: i) the stability analysis of these models is made much easier by mapping f to P and using the stability analysis of P, for which we know the exact expression and nature of its equilibria; ii) that expressing the parameters of P (α 1 , α 2 and α 3 ) as aggregates of the parameters of f helps to understand the effect of biological parameters on the systems equilibria -in particular through the graphical analysis of the α 2 -α 1 α 3 plane -and iii) that the models we consider all share a common structure despite the fact that they describe different mechanistic processes, which is encouraging for the transposability of studies made on one of these models.

S3.2.1 Noy-Meir model

The model we use in the main text of this study describes the dynamics of plants (whose biomass is written x) grazed by a constant population of consumer (whose maximal grazing rate is B), as coined by [START_REF] Noy-Meir | Stability of grazing systems: an application of predator-prey graphs[END_REF]. The variations of plant biomass are expressed by the following differential 36 equation:

dx dt := f (n) = rx 1 - x K - Bx A + x (5)
In the absence of grazing, the plant biomass grow logistically with a growth rate r up to a carrying capacity K (left term). On top of that, a fixed consumer harvest the plant population following a Holling type 2 harvesting rate (right term) where B is the maximal harvesting rate and

A is the biomass of plants at which the per plant biomass harvesting rate is half of its maximal value (half-saturation constant of the grazers): the grazing increases linearly with x at low plant biomass then reaches a plateau as each harvesting agent reaches its maximal harvesting rate.

Polynomial form: We rewrite f (x) as a quotient of polynomial:

f (x) = -r K x 3 + r(1 -A/K)x 2 + (Ar -B)x A + x = P(x) g(x)
Hence the coefficients of P(x) are:

Coefficient Expression α 1 Ar -B α 2 r(1 -A/K) α 3 -r/K
Table S3.1: Coefficient of the polynomial P(x) describing the behaviour of the Noy-Meir model.

The main control parameter -grazing intensity (B) -is only found in α 1 . Varying this parameter moves the system vertically in the α 2 -4α 1 α 3 plane (Fig. S3.2a): it goes from the domain with a single positive equilibrium (green area) to the domain with a single null equilibrium (extinction, red area) through the bistable domain (orange area) as B increases. Note that in the unlikely case where A > K (i.e., the grazing saturates at a biomass level greater than the carrying capacity), α 2 is negative so the system goes from a positive biomass to extinction without going through a bistable phase as B increases.

37

S3.2.2 Klausmeier model

This model coined by [START_REF] Klausmeier | Regular and irregular patterns in semiarid vegetation[END_REF] and describes catastrophic shifts in arid grassland, but here the bistability arises from positive feedbacks between plant biomass (x) and water availability (w):

dw dt = A -Lw -Rwx 2 dx dt = RJwx 2 -Mx
with the following paramters:

• A the water input

• L the water loss rate • R the rate at which plants uptake water

• J a conversion rate from water mass to plant biomass (i.e., how much plant biomass is created from 1 unit of water mass).

• M the plant mortality rate.

The water dynamics depend on a constant water input (A, e.g., the precipitation rate), water losses that are proportional to the amount of available water (Lw, e.g., evaporation and infiltration in deep soil layers) and a water uptake by plants (Rwx 2 ) that is proportional to both the amount of proportional to the amount of available water and the square of plant biomass (the square term representing the facilitation between plants allows the system to be bistable). The plant dynamics are described by a growth term that is proportional to the water uptake (RJwx 2 ) and a death rate proportional to plant biomass (Mx).

It is common to assume that water dynamics are way faster than plant dynamics in models describing 38 plant-water dynamics. We can thus calculate the water equilibrium w * by solving the equation dw dt = 0:

w * = A L + Rx 2
By replacing w by its equilibrium w * in the expression of water dynamics, we can describe the system using a single equation:

dx dt := f (x) = RJAx 2 L + Rx 2 -Mx
Polynomial form: Similarly to the Noy-Meir model, we can rewrite the equation above as a quotient of polynomials:

f (x) = -RMx 3 + ARJx 2 -LMx L + Rx 2 = P(x) g(x)
The coefficients of P(x) are:

Coefficient Expression α 1 -LM α 2 ARJ α 3 -RM
Table S3.2: Coefficient of the polynomial P(x) describing the behaviour of the Klausmeier model.

The two main control parameters here are the plant mortality (M) and the water input (A). The plant mortality is only present in the al pha 1 and al pha 3 so -similarly to the Noy-Meir model -varying M moves the system vertically in the α 2 -4α 1 α 3 plane (Fig. S3.2b): it goes from the frontier between the domain with a single positive equilibrium (green area) and the bistable domain (orange area) to the domain with a single null equilibrium (extinction, red area) as M increases.

The water input parameter (A) is only present in α 2 , so varying it moves the system horizontally in the α 2 -4α 1 α 3 plane (Fig. S3.2b): increasing A moves the system from the red domain (monostable, extinction) to the orange domain (bistable).

Note that for this model, α 1 and α 3 are always negative (hence 4α 1 α 3 is always positive) and α 2 is always positive. Thus the system is necessarily in the top-right quadrant of the α 2 -4α 1 α 3 plane and thus is either bistable or admits only the null equilibrium (it cannot reach the green domain in the α 2 -4α 1 α 3 plane).

S3.2.3 Rietkerk & Van de Koppel model

The third model of desertification that we examine was coined van de Koppel et al. (1997) and also describes the dynamics of soil water (W ) and plant biomass (x):

dw dt = I x + Aq x + A -Lw -Uwx dx dt = rx(1 -x/K) w w + B -Mx
Water (w) enters the soil through an infiltration facilitated by plants presence (I(x + As)/(x + A),

with I the water input) and exits it through abiotic losses (-LW ) and plant uptake (-Uwx). The plants (x) grow following a logistic term (rx(1x)/K) multiplied by a saturating function of w (w/(w + B)) describing how water affects plant growth; and they lose biomass at a constant rate M (-Mx). Similarly to the Klausmeier model, we can find the water equilibrium w * by solving dw/dt = 0:

w * = I(x + Aq) (L +Ux)(x + A)
We can inject this equilibrium (under the assumption that water dynamics are faster than plant dynamics) into the expression of dx/dt to find:

dx dt := f (x) = rx 1 - x K w * w * + B -Mx
Polynomial form: By developing this expression, we find that it is once again a quotient of polynomials:

40 f (x) = P(x) g(x)
with:

g(x) = A(Iq + BL) + x(I + B(L + AU)) + x 2 BU and:

P(x) = x 3 (-Ir/K -MBU) + x 2 (Ir(1 -Aq/K) -IM -MBL -MBUA) + x(IAq(r -M) -MBAL)
whose coefficients are:

Coefficient Expression α 1 IAq(r -M) -MBAL α 2 Ir(1 -Aq/K) -IM -MBL -MBUA α 3
-Ir/K -MBU Table S3.3: Coefficient of the polynomial P(x) describing the behaviour of the model described by Rietkerk and Van de Koppel.

This model is a bit more complicated than the previous ones and would warrant a lengthy analysis of all its parameters and how they affect the systems equilibria. However, we will here focus on two simple results:

Firstly, the polynomial form provides a simple way to analyse the system: while the model uses 9 parameters, we can reduce the complexity of its analysis through the three compound parameters α 1 , α 2 and α 3 . For any set of parameters, it is easy to compute the al pha and then i) to determine the number and nature of the equilibria through the α 2 -4α 1 α 3 plane (Fig. S3.2c) and ii) to determine the exact values of these equilibrium through the analytical expresions of r 1 and r 2 that we derived in the section S3.1.

Secondly, we can analyze the impact of varying a given parameter using the α 2 -4α 1 α 3 plane. For the set of parameters used in Fig. S3.2c, we can see that the two control parameters (mortality rate

S4 Sensitivity analysis on network connectivity

We conducted a sensitivity analysis on network conductivity by conducting the same simulations as in the main text but on networks with higher or lower connectivity (see S2 for the generating methods). We generated regular graphs where each patch had exactly 2 or 8 neighbours, which respectively results in 100 patches on a circle with connections to the 2 nearest neighbours (1-D system, mean degree = 2) and a 10*10 lattice with periodic boundary conditions and connections to the 8 nearest neighbours (2-D system, mean degree = 8). We generated Erdös-Rényi graphs with a target connectivity of 2 or 8, but because we constrained them to be connected graphs, this resulted in an average connectivity slightly higher than 2 (2.7 on average) for the low connectivity networks.

We also generated RGGs with a target connectivity of 8 (which resulted in a realized connectivity of 7.9 on average), but couldn't make RGGs with a lower connectivity than in the main text because of constraints on the generating process.

S4.1 General patterns

The general patterns observed in the main text stayed qualitatively unchanged (Fig. S4.1 and S4.2):

hysteresis size first decrease with dispersal rate (until µ ∼ 0.1) and then increased to a plateau.

The effect of landscape type on hysteresis was conserved for highly connected graphs (regular graphs > Erdös-Rényi graphs > RGGs, fig. S4.2) We suspect that this is due to differences in degree distributions: in regular graphs, all patches have the same degree (4 or 8 neighbours), which explains their larger hysteresis. Erdös-Rényi networks and RGGs have some variance in their node degrees, so the perturbations were likely to happen on a low degree node and to initiate a large scale shift, explaining their hysteresis size in between that of regular graphs and of OCNs. The difference between RGGs and Erdös-Rényi graphs could be also due to differences degree distribution: while they have a similar average degree, RGGs have a wider degree distribution, making it more likely that a pertubation happens on a low-degree node than in Erdös-Rényi graphs, hence their smaller hysteresis.

S5 Analysis of a two-patch system

At high dispersal (Fig. S5.1 a, b & c), the densities at equilibrium are the same in the two patches because of the strong homogenization between patches: all equilibria are situated on the x 1 = x 2 plane. Spatial bistability -the stable coexistence of a high-and low-biomass states in the two patches -is not possible, and biomass differences between the two patches can only exist as a transient state towards an homogeneous stable state. The whole system behaves as a single bistable unit with two stable equilibria ((x 1 , x 2 ) = (0, 0) and (x 1 , x 2 ) = (r 1 , r 1 )) separated by an unstable equilibrium ((x 1 , x 2 ) = (r 2 , r 2 )) over the same range of harvesting rate as for an isolated patch (50 < B < 56.25).

As dispersal decreases, the equilibria on the x 1 = x 2 plane keep their values and nature, but branches appear outside of the x 1 = x 2 plane at high harvesting rates (right before the tipping point toward the desert equilibrium) (Fig. S5.1 d, e & f). These branches appear from the unstable branch ((x 1 , x 2 ) = (r 2 , r 2 )) that separates the desert equilibrium ((x 1 , x 2 ) = (0, 0)) from the vegetated equilibrium ((x 1 , x 2 ) = (r 1 , r 1 )) and are themselves unstable. At really low dispersal, spatial bistability becomes possible (meaning the coexistence of patches with different plant biomasses in the metacommunity) since part of the branches outside the x 1 = x 2 plane become stable (Fig. S5.1 g, h & i). These new stable equilibria are associated with strong differences in biomass between the two patches, with one patch at a biomass really close to r 1 while the other is close to 0.

F1Figure 1 :

 1 Figure 1: Bifurcation diagrams of a single patch (a) and of two patch systems (b, c) at high (b, µ = 0.3) and low (c, µ = 0.005) dispersal rates. Full lines are stable equilibria, dashed lines are unstable equilibria. (a) Equilibria of the biomass (x) of an isolated patch as a function of the grazing pressure (B). (b, c) Equilibria of the biomasses (x 1 , x 2 ) of two patches connected by dispersal, as a function of the grazing pressure (B).

Figure 2 :

 2 Figure 2: State diagram of a regular grid made of 10*10 bistable patches as a function of the harvesting rate (B) and the dispersal rate (mu) (a, b, c) and bifurcation diagrams of the average biomass for various dispersal rates (d, e, f). (a) and (b) state diagrams of a regular grid of 10*10 bistable patches when starting from 5 random occupied patches in a otherwise empty landscape (a) or 5 random empty patches in an otherwise occupied landscape (b). The colors denote the final state of the system: "fully occupied" (green) means that the total equilibrium biomass was greater than 99% of the maximal equilibrium biomass, "empty" means that the total equilibrium biomass was smaller than 1% of the maximal equilibrium biomass, and "partially occupied" means that the system was neither empty nor fully occupied. (c) general state diagram of the aforementioned system, obtained by superimposing panels (a) and (b). (d), (e) and (f) bifurcation diagrams of the average biomass of a regular grid (10*10 bistable patches) for different dispersal rates ((d): µ = 0.5, (e): µ = 0.05, (f): µ = 0.001). State space limits from panels (a), (b) and (c) were obtained as the average of 50 replicates; panels (d), (e) and (f) are the bifurcation diagram of a single replicate each.

Figure 3 :

 3 Figure 2: State diagram of a regular grid made of 10*10 bistable patches as a function of the harvesting rate (B) and the dispersal rate (mu) (a, b, c) and bifurcation diagrams of the average biomass for various dispersal rates (d, e, f). (a) and (b) state diagrams of a regular grid of 10*10 bistable patches when starting from 5 random occupied patches in a otherwise empty landscape (a) or 5 random empty patches in an otherwise occupied landscape (b). The colors denote the final state of the system: "fully occupied" (green) means that the total equilibrium biomass was greater than 99% of the maximal equilibrium biomass, "empty" means that the total equilibrium biomass was smaller than 1% of the maximal equilibrium biomass, and "partially occupied" means that the system was neither empty nor fully occupied. (c) general state diagram of the aforementioned system, obtained by superimposing panels (a) and (b). (d), (e) and (f) bifurcation diagrams of the average biomass of a regular grid (10*10 bistable patches) for different dispersal rates ((d): µ = 0.5, (e): µ = 0.05, (f): µ = 0.001). State space limits from panels (a), (b) and (c) were obtained as the average of 50 replicates; panels (d), (e) and (f) are the bifurcation diagram of a single replicate each.

Figure 4 :

 4 Figure4: Effect of different landscape types (spatial structures) on the regional dynamics of metapopulations (purple: regular graphs; green: Erdös-Rényi graphs; brown: random geometric graphs (RGGs); blue: optimal channel network (OCNs)). (a) Size of the hysteresis (computed as the area between the upper and lower branches of the bifurcation diagram of the average biomass) as a function of the dispersal rate. (b) and (c) large scale tipping points as a function of the dispersal rate: value of the harvesting rate (B) at which the system loses (b) or gains (c) the most biomass. All quantities are computed as the mean (full line) and standard deviation (colored areas) of 50 replicates.

Figure 5 :

 5 Figure5: Effects of different perturbation types on the regional dynamics of metapopulations (full line: clustered perturbations, dashed lines: dispersed perturbations) and different landscape types (purple: regular grid; brown: random geometric graphs (RGGs); blue: optimal channel networks (OCNs)). (a) Size of the hysteresis as a function of the dispersal rate. (b) and (c) regional scale tipping points as a function of the dispersal rate: value of the harvesting rate (B) at which the system loses (b) or gains (c) the most biomass. All quantities are computed as the mean (lines) and standard deviation (colored areas) of 50 replicates.

Figure 6 :

 6 Figure6: Effect of different perturbation types on the regional dynamics of the metapopulation (full line: perturbation of high-degree patches, dashed lines: perturbation of low-degree patches) and different landscape types (green: Erdös-Rényi graphs; brown: random geometric graphs (RGG); blue: optimal channel networks (OCN)). (a) Size of the hysteresis as a function of the dispersal rate. (b) and (c) regional scale tipping points as a function of the dispersal rate: value of the harvesting rate (B) at which the system loses (b) or gains (c) the most biomass. All quantities are computed as the mean (lines) and standard deviation (colored areas) of 50 replicates.

Figure

  Figure S3.1: Description of P behaviour. (a) Phase diagram of P depending on α 2 and 4α 3 α 1 . In the red area, 0 is the only equilibrium in R + and the system always goes to extinction (see panel c). In the orange area, the system has three equilibrium in R + ([0, r 1 , r 2 ]) and goes either to extinction or r 2 depending on the initial x (see panel d). In the green area, 0 is an unstable equilibrium and the system always goes to r 1 that is stable (see panel e). The full line shows a transect happening when increasing both α 2 from 0.8 to 2.4 and α 1 from -0.6 to 0.2 (α 3 is fixed at -0.5) and the dots show where the panels c, c and e where taken. (b) bifurcation diagram (equilibria of P(x)) along the transect represented in the panel a (the x axis figures the position along the line: 0 is the point (c) and 1 is the point (e)). (c) a monostable parameterization of P(x) resulting in extinction (parameters values of the point (c) in panel a: α 1 = -0.6; α 2 = 0.8; α 3 = -0.5. (d) a bistable parameterization of P(x) (parameters values of the point (d) in panel a: α 1 = -0.4; α 2 = 1.2; α 3 = -0.5. (e) a monostable parameterization of P(x) resulting in a positive equilibrium (parameters values of the point (e) in panel a: α 1 = 0.2; α 2 = 2.4; α 3 = -0.5.

Figure

  Figure S3.1: Description of P behaviour. (a) Phase diagram of P depending on α 2 and 4α 3 α 1 . In the red area, 0 is the only equilibrium in R + and the system always goes to extinction (see panel c). In the orange area, the system has three equilibrium in R + ([0, r 1 , r 2 ]) and goes either to extinction or r 2 depending on the initial x (see panel d). In the green area, 0 is an unstable equilibrium and the system always goes to r 1 that is stable (see panel e). The full line shows a transect happening when increasing both α 2 from 0.8 to 2.4 and α 1 from -0.6 to 0.2 (α 3 is fixed at -0.5) and the dots show where the panels c, c and e where taken. (b) bifurcation diagram (equilibria of P(x)) along the transect represented in the panel a (the x axis figures the position along the line: 0 is the point (c) and 1 is the point (e)). (c) a monostable parameterization of P(x) resulting in extinction (parameters values of the point (c) in panel a: α 1 = -0.6; α 2 = 0.8; α 3 = -0.5. (d) a bistable parameterization of P(x) (parameters values of the point (d) in panel a: α 1 = -0.4; α 2 = 1.2; α 3 = -0.5. (e) a monostable parameterization of P(x) resulting in a positive equilibrium (parameters values of the point (e) in panel a: α 1 = 0.2; α 2 = 2.4; α 3 = -0.5.

  Figure S4.2: Sensitivity analysis: networks with an average degree of 8. Characterization of the large scale dynamics for different spatial structures of metapopulation (purple: regular grid; green: Erdös-Rényi graphs; brown: random geometric graphs). (a) Size of the hysteresis (computed as the area between the upper and lower branches of the bifurcation diagram of the average biomass) as a function of the dispersal rate. (b) and (c) large scale tipping points as a function of the dispersal rate: value of the harvesting rate (B) at which the system loses (b) or gains (c) the most biomass. All quantities are computed as the mean (full line) and standard deviation of 50 replicates.

  Figure S6.1: Examples of three landscapes used in the main text analysis. (a) An Erdös-Rényi graph, (b) a random geometric graph (RGG) and (c) an optimal chanel network (OCN).

Figure S6. 2 :b

 2 Figure S6.2: Examples of a low-connectivity Erdös-Rényi graph used in the sensitivity analysis.

  for which the dispersal is strong enough for a local restoration to affect neighbouring patches, but still sufficiently weak to avoid the dilution of the added biomass in space, triggering a travelling wave of restoration. Interestingly, we observed a large hysteresis size at both very low and high dispersal rates, but with opposite underlying mechanisms: in the case of low dispersal, large scale shifts are hard to induce because local perturbations in a few patches won't affect the neighbouring patches at all. At high dispersal, large scale shifts are also hard to induce, but because local perturbations quickly get diluted across the whole landscape.It is worth noting that these variations with dispersal rate are reminiscent of[START_REF] Zelnik | The three regimes of spatial recovery[END_REF]'s regimes of spatial recovery: local restoration efforts could only spread in space at intermediate dispersal rates (what Zelnik called "rescue recovery"), but not at very low dispersal ("isolated recovery": local restoration does not affect neighbouring patches) or at very high dispersal ("mixing recovery": the local restoration effort is diluted over the whole landscape).

Supplementary material S1 Metapopulation model

We used a model of spatially structured population describing the dynamics of n patches linked by the dispersal of individuals. The dynamics of the biomass (x i ) inside the patch i are described by the following differential equation:

where f (x i ) is a function describing the local dynamics of the patch and g(x i , x • ) describes the net dispersal between the patch i and its neighbors.

Local dynamics We describe the local dynamics with a model coined by [START_REF] Noy-Meir | Stability of grazing systems: an application of predator-prey graphs[END_REF] that exhibits bistability and catastrophic shifts between a null-biomass state (x i = 0) and a positive equilibrium (x i = x * ) (Eq. 3). It describes the dynamics of a biotic resource x i (e.g., the biomass of plants) that grows logistically with a growth rate r and a carrying capacity K. This resource is harvested/consumed by a fixed consumer with a Holling type II functional response (with B the maximal harvesting rate and A the resource biomass at which the harvesting is half of B):

Although it was first coined to describe the dynamics of a plant biomass is a pasture grazed by herbivores [START_REF] Noy-Meir | Stability of grazing systems: an application of predator-prey graphs[END_REF], this model is fairly general as logistic growth is commonly used to model the growth of a wide range of organisms (but see [START_REF] Mallet | The struggle for existence. how the notion of carrying capacity, k, obscures the links between demography, darwinian evolution and speciation[END_REF]), and type II functional responses are the norm for numerous consumers [START_REF] Jeschke | Consumer-food systems: why type i functional responses are exclusive to filter feeders[END_REF][START_REF] Kalinkat | Body masses, functional responses and predator-prey stability[END_REF]. This system admits two tipping points (B = AR and B = Ar + K(1 -AK)/4). It is bistable in this range (AR < B < Ar + K(1 -AK)/4) and admits two stable equilibria (x i = 0 and x i = r 1 ) separated by an unstable equilibrium (0 < r 2 < r 1 ).

M and water input I) have a similar effect on the system. The system moves from a stable positive equilibrium (green area) towards extinction (red area) through a bistable state (orange area) when increasing the mortality (M, red line in Fig. S3.2c) or decreasing the water input (I, blue line).

S3.2.4 Logistic growth with an allee effect

Lastly, we focus on the model of a logistic growth with an Allee effect. For simplicity, we use a model where x represents the species density rescaled to the carrying capacity (i.e., x is the species density in proportion of the carrying capacity). This system is described by the following equation:

where r is the growth rate and c is the Allee effect threshold, i.e., the density under which the species growth is negative.

Polynomial form: Here, the equation describing the system is already a third degree polynomial (e.g., g(x) = 1) that we can express by developing the expression above:

Hence the coefficients of the polynomial form are:

Table S3.4: Coefficient of the polynomial P(x) describing a logistic growth with an Allee effect.

Because this model is very simple by design, the polynomial analysis is not very useful here. It is however interesting to note that it fits into the same framework as the more mechanistic models presented above. The figures S3.2d shows how this model is situated in the α 2 -4α 1 α 3 plane as c varies.

Hysteresis size was similar accross landscapes types for lowly connected graphs (Fig. S4.1), which we conjecture is due to the very low degree heterogeneity at low connectivity.

Lastly, hysteresis size was smaller in less connected networks (Fig. S4.1) and higher in more connected networks (Fig. S4.1), which is consistent with our interpretation that biomass dilution in highly connected networks prevent the spread of local shifts.

S4.2 Perturbation modality

Lastly, the effect of perturbation modality was similar to what we observe in the main text. Targeting clustered patches resulted in a smaller hysteresis than targeting patches dispersed over the landscape in both lowly and highly connected networks (Fig. S4.3 and S4.4).

Targeting high or low degree node had only a marginal effect in lowly-connected Erdös-Rényi networks and no effect otherwise (Fig. S4.5 and S4.6), similarly to what we observed in 4-degree networks.


Maryse Vanderplanck
email: maryse.vanderplanck@cefe.cnrs.fr

Benoit Lapeyre
email: benoit.lapeyre@cefe.cnrs.fr

Shéhérazade Lucas
email: sheherazade.lucas@outlook.com

Magali Proffit
email: magali.proffit@cefe.cnrs.fr

Francesco Parisi

Klaus H Hoffmann

Enrico Ruzzier

Simone Sabatelli

Ozone Induces Distress Behaviors in Fig Wasps with a Reduced Chance of Recovery

Keywords: Agaonidae, air pollution, abiotic risk, pollinator, behavioral change, recovery

Simple Summary: Ecological interactions among organisms underpin the stability of ecological networks, which are responsible for species biodiversity in ecosystems. These interactions are currently threatened by environmental risks, mainly due to human activities, such as air pollution. Among air pollutants, tropospheric ozone (O 3 ) is known to disrupt chemical communication between plants and their pollinators. Alarmingly, its concentration is likely to increase by two-four-fold in the next two decades. However, the direct effects of O 3 on the behavior of pollinators themselves have not been investigated so far, even though insect behavior is key to their ecological interactions. In this study, we evaluated the potential effects of O 3 at different field-realistic concentrations on the behavior of the fig wasp Blastophaga psenes, the exclusive pollinator of the Mediterranean fig species Ficus carica. We found that O 3 , even at low concentrations, induced abnormal motility in fig wasps, and that exposed individuals might only have a reduced chance of recovery. Overall, our findings indicate that O 3 can affect pollinator behavior, which may have detrimental implications for pollination systems.

Introduction

Current global change substantially threatens the balance of ecosystems through decreased species abundance and diversity [START_REF] Potts | Global Pollinator Declines: Trends, Impacts and Drivers[END_REF]. Many synergistic drivers contribute to this loss in biodiversity, including habitat destruction/fragmentation [START_REF] Goulson | Decline and Conservation of Bumble Bees[END_REF][START_REF] Winfree | A Meta-Analysis of Bees'Responses to Anthropogenic Disturbance[END_REF], the use of agrochemicals [START_REF] Kevan | Log-Normality of Biodiversity and Abundance in Diagnosis and Measuring of Ecosystemic Health: Pesticide Stress on Pollinators on Blueberry Heaths[END_REF][START_REF] Halm | Modes of Honeybees Exposure to Systemic Insecticides: Estimated amounts of Contaminated Pollen and Nectar Consumed by Different Categories of Bees[END_REF], decreasing resource diversity [START_REF] Biesmeijer | Parallel Declines in Pollinators and Insect-Pollinated Plants in Britain and the Netherlands[END_REF], climate change [START_REF] Williams | Can Vulnerability Among British Bumblebee (Bombus) Species Be Explained by Niche Position and Breadth?[END_REF][START_REF] Dormann | Prediction Uncertainty of Environmental Change Effects on Temperate European Biodiversity[END_REF], and air pollution [START_REF] Agathokleous | Ozone Affects Plant, Insect, and Soil Microbial Communities: A Threat to Terrestrial Ecosystems and Biodiversity[END_REF][START_REF] Blande | Effects of Air Pollution on Plant-Insect Interactions Mediated by Olfactory and Visual Cues[END_REF].

All these environmental risks impact living organisms through changes in their phenology, distribution, and behavior [START_REF] Wong | Behavioral Responses to Changing Environments[END_REF][START_REF] Gérard | Global Warming and Plant-Pollinator Mismatches[END_REF]. Such processes depend on species' ecological traits and often involve cascading effects [START_REF] Furlan | Small Population Size and Extremely Low Levels of Genetic Diversity in Island Populations of the Platypus, Ornithorhynchus anatinus[END_REF]. Among these effects, changes in behavior may lead to drastic modifications and/or disruptions of crucial species interactions, including plant-pollinator interactions [START_REF] Wong | Behavioral Responses to Changing Environments[END_REF].

Among pollinators, insects exhibit spectacular and diverse behaviors that are key to their participation in ecological interactions. These behaviors can be acutely sensitive to environmental stressors, such as exposure to chemicals [START_REF] Augusiak | The Influence of Insecticide Exposure and Environmental Stimuli on the Movement Behaviour and Dispersal of a Freshwater Isopod[END_REF]. Specifically, a change in behavior is very often the first response of an organism to an environmental change. Such behavioral modifications can potentially facilitate adaptation in a changing world by improving the survival and reproduction of species, or hinder it by interfering with physiological processes or disrupting species interactions (reviewed in [START_REF] Wong | Behavioral Responses to Changing Environments[END_REF]). These interactions between species underpin the stability of the ecological networks that shape biodiversity. Understanding changes in the behavior of interacting species is therefore essential to predict the future stability of ecological networks in a changing world. However, the effects of environmental stressors on pollinator behaviors remain understudied, especially the effect of air pollution. Tropospheric ozone (O 3 ) is one of the most damaging air pollutants, known to have detrimental effects on organisms [START_REF] Iriti | Oxidative Stress, the Paradigm of Ozone Toxicity in Plants and Animals[END_REF][START_REF] Vanderplanck | Ozone Pollution Alters Olfaction and Behavior of Pollinators[END_REF], and whose concentration is likely to increase by two-four-fold in the next two decades [START_REF] Cooper | Global Distribution and Trends of Tropospheric Ozone: An Observation-Based Review[END_REF]. Here, we tested whether ozone exposure episodes (local seasonal peaks of >40 ppb

Materials and Methods

Biological Model

Blastophaga psenes (Hymenoptera, Agaonidae) is a solitary and tiny fig wasp, which is diurnal (active during the day), lives generally less than one day, and does not feed at the adult stage [START_REF] Kjellberg | Longevity of a Fig Wasp (Blastophaga psenes)[END_REF]. This fig wasp species is intimately associated with Ficus carica (Moraceae) for its reproduction and is its exclusive pollinator. Blatophaga psenes has two generations per year that coincide with flowering of fig trees. Wasp larvae overwinter as the last larval stage in the male figs produced in summer. In spring, the overwintering wasps complete their development and emerge in summer in search of receptive figs (usually between 10 am and 5 pm), generally located on a different tree, in which they multiply [START_REF] Kjellberg | Longevity of a Fig Wasp (Blastophaga psenes)[END_REF][START_REF] Proffit | Chemical Signal is in the Blend: Bases of Plant-Pollinator Encounter in a Highly Specialized Interaction[END_REF]. Because of their ecological traits, fig wasps are probably highly susceptible to environmental stressors, as they have a reduced chance of recovery (i.e., reduced lifespan) and plasticity (i.e., reduced lifespan and highly host-specialized interaction), and no possibility for nutritional resilience (i.e., an inability to feed as adults).

Ozone Exposure and Behavioral Assays

In July 2018, newly emerging adult female wasps were collected from mature figs, taken randomly from different male trees ("Terrain d'Expérience" of the "Centre d'Écologie Fonctionnelle et Évolutive-CEFE", pesticide-free, Montpellier, France), and tested shortly after their exit from their natal fig. Groups of fig wasps (around 11 individuals) were placed into a laboratory fumigation chamber, which consisted of a glass bottle of 500 mL with a filter paper of 2 cm × 2 cm loaded with 200 µL of distilled water prior to O 3 exposure. Ozone was produced using the photolysis of molecular oxygen subjected to UV radiation at a wavelength of 185 nm (UV photometric Ozone Analyzer, Model 49i, Thermo Scientific) and delivered continuously in a flow through the fumigation chamber (flow rate of 1.5 L•min -1 ). One extremity of the fumigation chamber was connected to an analyzer generator in the generator mode, pushing air containing different O 3 concentrations into the bottle at a flow rate of 1.5 L•min -1 , while the other extremity was connected to an analyzer generator in the analyzer mode, where air was extracted at a flow rate of 1.5 L•min -1 to ensure that the desired O 3 concentration was present in the bottle. Before entering the system, the air was cleaned of any VOCs by using an activated carbon filter. We used exclusively Teflon tubes to connect the pump, the VOC filter, and the O 3 generator and analyzer (Figure 1

Statistical Analyses

All analyses were performed in R version 3.4.0 [START_REF] Team | A Language and Environment for Statistical Computing[END_REF]. To estimate the probability of an individual displaying motility after an O 3 episode, non-parametric Kaplan-Meier curves were estimated using the "survfit" function (R package "survival"). Individuals that displayed abnormal motility or distress behaviors (including death) were the exact observations, while those that remained active or were not followed-up until the end of the survival tracking were the censored data. We used a log-rank test ("survdiff" function, R package "survival") to test for differences in Kaplan-Meier curves among the O 3 concentrations. Post hoc tests were conducted using pairwise log-rank tests adjusted for false discovery rate (FDR) ("pairwise_survdiff" function, R package "survminer"). Furthermore, we used multivariate generalized linear models (ManyGLM) to analyze the effects of O 3 exposure on fig wasp behavior in detail ("mvabund" R package [START_REF] Wang | Mvabund Statistical Methods for Analysing Multivariate Abundance Data[END_REF]). By correcting for negative binomial responses, this multivariate method allows for the handling of zero-inflated multiple count data without normality assumption [START_REF] Wang | Mvabund-An R Package for Model-Based Analysis of Multivariate Abundance Data[END_REF][START_REF] Warton | Distance-Based Multivariate Analyses Confound Location and Dispersion Effects[END_REF]. Moreover, it has more statistical power than distance-based multivariate analyses such as perMANOVA [START_REF] Warton | Distance-Based Multivariate Analyses Confound Location and Dispersion Effects[END_REF][START_REF] Warton | Regularized Sandwich Estimators for Analysis of High-Dimensional Data Using Generalized Estimating Equations[END_REF]. As there are no extant mixed-effect models for the ManyGLM method, we analyzed the effect of O 3 exposure for each hour of survival tracking separately. When significant difference was detected, multiple pairwise comparisons were conducted followed by the univariate test procedure implemented in ManyGLM. We reported unadjusted p-values for all tests, because the detectability after adjustment for multiple comparisons would be very low given the high number of comparisons, but interpreted the results with caution [START_REF] García | Escaping the Bonferroni Iron Claw in Ecological Studies[END_REF][START_REF] Nakagawa | A Farewell to Bonferroni: The Problems of Low Statistical Power and Publication Bias[END_REF].

Results and Discussion

The probability of an individual displaying motility after an O Such negative effects on fig wasps' behavior could be explained as due to oxidative stress triggered by O 3 episodes. Previously, O 3 has been shown to induce molecular damage such as the oxidation of proteins, lipid peroxidation, and damage to DNA; it also causes the deregulation of intracellular signal transduction, which could disrupt the whole organism and lead to death (reviewed in [START_REF] Kodrík | Hormonal Regulation of Response to Oxidative Stress in Insects-An Update[END_REF]). To support this hypothesis, further experiments are required to evaluate the direct effect of O 3 on fig wasp individuals by assessing oxidative stress biomarkers, such as levels of lipid peroxidation and total reactive antioxidant potential [START_REF] Kodrík | Hormonal Regulation of Response to Oxidative Stress in Insects-An Update[END_REF][START_REF] Telesnicki | Direct Effect of Ozone Pollution on Aphids: Revisiting the Evidence at Individual and Population Scales[END_REF]. However, as fig wasps probably vary in their susceptibility to O 3 (i.e., among-individual variation), differences in these oxidative stress markers should be assessed in several individuals every hour during the hourly monitoring, which involves an experimental setup with adequate replicates. Indeed, comparing the effect of exposure regimes by assessing oxidative stress markers only in surviving individuals at the end of the trials (i.e., after five hours) would lead to biases, such as the absence of difference, as only the tolerant fig wasps would be accounted for (i.e., oxidative stress underestimated) [START_REF] Telesnicki | Direct Effect of Ozone Pollution on Aphids: Revisiting the Evidence at Individual and Population Scales[END_REF]. With respect to the chances of recovery of O 3 -exposed fig wasps, our study noted that the effect of exposure lasted up to five hours after the cessation of O 3 exposure, at least, with even stronger effects over time. This limited possibility of recovery is reinforced by the ecological traits of this pollinator species. Indeed, the flight period of fig wasps (i.e., diurnal species, active between 10 a.m. and 5 p.m.) exposes them to high O 3 concentrations. As an additional point, their reduced lifespan (generally less than 24 h), as well as their inability to feed as adults, does not allow for the activation of endogenous antioxidant defense systems (i.e., energetically costly processes requiring time to be activated) [START_REF] Beaulieu | Behavioural Antioxidant Strategies to Cope with High Temperatures: A Study in a Tropical Butterfly[END_REF] and prevents them from benefiting from a protective effect of dietary antioxidants [START_REF] Bonvehí | Evaluation of Polyphenolic and Flavonoid Compounds in Honeybee-Collected Pollen Produced in Spain[END_REF][START_REF] Johnson | Plant Phenolics as Dietary Antioxidants for Herbivorous Insects: A Test with Genetically Modified Tobacco[END_REF].

By Motility 11 ± 3 11 ± 3 10 ± 3 10 ± 3 9 ± 3 9 ± 3 Abnormal motility 0 ± 0 0 ± 0 1 ± 2 1 ± 2 0 ± 0 0 ± 0 Moderate distress behavior 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 Severe distress behavior 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 Death 0 ± 0 0 ± 0 0 ± 1 0 ± 1 0 ± 1 0 ±

± 0 0 ± 0 1 ± 1 1 ± 1 2 ± 2 3
± 2 Moderate distress behavior 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 Severe distress behavior 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 Death 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0

± 0 0 ± 0 1 ± 1 1 ± 1 2 ± 2 2 ± 3 Moderate distress behavior 0 ± 0 0 ± 0 1 ± 2 1 ± 1 1 ± 1 1 ± 2 Severe distress behavior 0 ± 0 0 ± 0 0 ± 1 1 ± 3 0 ± 0 1 ± 2 Death 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 1 1 ± 1

  O 3 ) result in a change in the behavior of the fig wasp Blastophaga psenes, the exclusive pollinator of the Mediterranean fig species Ficus carica, and discuss the implications for the fig/fig wasp pollination system.

  ). Individuals were exposed to 0 (control, n = 10 with 11 ± 3 individuals per trial), 40 (n = 10 with 11 ± 2 individuals per trial), 80 (n = 11 with 11 ± 3 individuals per trial), 120 (n = 9 with 12 ± 1 individuals per trial), or 200 ppb (n = 10 with 11 ± 3 individuals per trial) for 60 min based on the recorded O 3 episodes in the French Mediterranean region [17,20,21]. All O 3 exposures were conducted in a greenhouse ("Terrain d'Expérience") maintained at 25 • C under natural light, generally between 11 a.m. and 12 a.m.

Figure 1 .

 1 Figure 1. Experimental setup for the controlled and continuous flow during O 3 exposure of fig wasps at 0, 40, 80, 120, or 200 ppb. After O 3 exposure, the fumigation chambers were disconnected from the O 3 generator and analyzer and maintained at 25 • C under natural light inside the greenhouse. We then performed observations to assess the probability of motility of the exposed fig wasps, without disturbing them inside the glass bottle. Groups of individuals were monitored hourly for 5 h to track survival and record behavioral responses. Five distinct behaviors were observed: (1) motility, (2) abnormal motility, (3) moderate distress behavior (i.e., individuals alternately lying on their backs or upright), (4) severe distress behavior (i.e., individuals lying on their backs, unable to right themselves), and (5) death.

  3 episode significantly differed according to the O 3 concentration (Kaplan-Meier analysis, χ 2 = 180, df = 4, p < 0.001). While control individuals (exposed to 0 ppb O 3 ) remained mostly active over time and even until the end of the tracking (after five hours), the individuals exposed to O 3 exhibited a decreasing probability of being active over time, which was significant at 80, 120, and 200 ppb O 3 (p < 0.001 for 0 vs. 80 ppb O 3 , 0 vs. 120 ppb O 3 and 0 vs. 200 ppb O 3 ) with an even more pronounced decrease at 120 and 200 ppb O 3 compared to 80 ppb O 3 (p < 0.001 for 80 vs. 120 ppb O 3 and 80 vs. 200 ppb O 3 ) (Figure 2, Appendix A).

Figure 2 .

 2 Figure 2. Kaplan-Meier curves depicting the probability of motility over time for each O 3 regime. Crosses indicate points of censored data. Different letters indicate significant differences (p < 0.05) in the probability of motility between O 3 regimes based on pairwise log-rank tests. The multivariate analyses showed that the proportion of individuals displaying abnormal motility significantly increased one hour after O 3 exposure at 120 and 200 ppb compared to both the control (0 ppb) and 40 ppb (0 vs. 120 ppb O 3 , p = 0.011; 40 vs. 120 ppb O 3 , p = 0.016; 0 vs. 200 ppb O 3 , p = 0.042; and 40 vs. 200 ppb O 3 , p = 0.021) (Figure 3b, Appendix A). At three hours after the O 3 episode the proportion of individuals displaying motility significantly decreased when exposed to 200 ppb O 3 compared to the

Figure 3 .

 3 Figure 3. Behavior of fig wasps exposed to O 3 regimes for one hour (0, 40, 80, 120, and 200 ppb). After O 3 exposure, individuals were left in climate-controlled conditions and monitored hourly for 5 h to track survival and record behavioral responses. Each subfigure (a-f) corresponds to one of the monitoring moments. Five distinct behaviors were observed: (1) motility, (2) abnormal motility, (3) moderate distress behavior, (4) severe distress behavior, and (5) death. Data are expressed as the mean percentage of individuals displaying the different behaviors. Asterisks indicate significant differences (**, very significant, *** extremely significant) in the behavior of exposed fig wasps between O 3 regimes.

1 Exposure at 40 ppb O 3 (n = 10 with 11 ± 2 Individuals per Trial)

		0 h	1 h	2 h	3 h	4 h	5 h
	Motility	11 ± 2	11 ± 2	10 ± 2	10 ± 2	9 ± 3	8 ± 3
	Abnormal motility	0

Exposure at 80 ppb O 3 (n = 11 with 11 ± 3 Individuals per Trial)

		0 h	1 h	2 h	3 h	4 h	5 h
	Motility	11 ± 3	11 ± 3	9 ± 3	8 ± 3	8 ± 2	7 ± 3
	Abnormal motility	0

Acknowledgments:

We thank Candice Dubuisson for her valuable help during the experiment, as well as the personnel of the Platform TE (LabEx CeMEB, the French National Research Agency (ANR) under the "Investissements d'avenir"). We would like to thank T.J. Wood for proofreading and improvement of the English. We are further grateful to the reviewers whose comments improved the presentation of our study.

Funding: This research was partly supported by the programs International Research Project (IRP)-CNRS-MOST, ANR under the "Investissements d'avenir" program (reference ANR-16-IDEX-0006) and the French National Research Program for Environmental and Occupational Health of ANSES (reference 2018/1/138). Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable.

Data Availability Statement: Data available on request from the authors.

Contributions: Conceptualization, M.P. and B.L.; data acquisition, B.L. and S.L.; data analysis, M.V.; visualization, M.V.; writing-original draft preparation, M.V.; writing-review and editing, M.V. and M.P.; supervision, M.P.; funding acquisition, M.P. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.