Pierre Quévreux

Bart Haegeman

Michel Loreau

Spatial heterogeneity of interaction strength has contrasting effects on synchrony and stability in trophic metacommunities

Keywords: source-sink, stochastic perturbations, food chain, dispersal, asymmetry, conservation

Introduction

Since [START_REF] May | Will a large complex system be stable?[END_REF] demonstrated that stability was not an inherent property of ecological interaction networks, ecologists have been relentlessly looking for the mechanisms ensuring ecosystem stability. Spatial heterogeneity has long been identified as one of the main factors promoting the mechanisms underlying the maintenance of biodiversity and the stability of ecosystems. For instance, in competitive metacommunity models, spatial heterogeneity provides local favourable conditions to each species of the regional pool [START_REF] Holt | Spatial heterogeneity, indirect interactions, and the coexistence of prey species[END_REF][START_REF] Chesson | Mechanisms of maintenance of species diversity[END_REF][START_REF] Amarasekare | Spatial heterogeneity, source-sink dynamics, and the local coexistence of competing species[END_REF], which in turn ensures species persistence in less favourable patches by source-sink dynamics (Mouquet andLoreau, 2002, 2003;[START_REF] Loreau | Biodiversity as spatial insurance in heterogeneous landscapes[END_REF]. The stability of the temporal dynamics of species biomass is ensured by the asynchrony of the dynamics between patches, which leads to compensatory dycorresponding author: pierre.quevreux@cri-paris.org namics [START_REF] Loreau | Biodiversity as spatial insurance in heterogeneous landscapes[END_REF][START_REF] Loreau | Species synchrony and its drivers: Neutral and nonneutral community dynamics in fluctuating environments[END_REF]. In trophic metacommunities, spatial heterogeneity has also been identified as a stabilising factor [START_REF] Steele | Spatial heterogeneity and population stability[END_REF][START_REF] Hastings | Spatial heterogeneity and the stability of predator-prey systems[END_REF][START_REF] Hastings | Spatial heterogeneity and the stability of predator-prey systems: Predator-mediated coexistence[END_REF], but the underlying mechanisms are more complex due to the interplay between trophic and spatial dynamics.

Inspired by the description of fast and slow energy channels by soil ecologists (i.e., in terms of biomass turnover), Rooney et al. (2006) noted the stabilising effect of the asymmetry of energy flows in ecosystems with a food web model consisting of one mobile predator feeding on two energy channels. In their model, the asymmetry of energy flow is generated by different interaction strengths between predators and prey (i.e. increased attack rate in one energy channel compared to the other one, see Figure 1) and different consumption rates of a common resource by the two basal species, which in turn promotes the asynchrony of prey biomass dynamics in response to perturbations. Although synchrony patterns are tightly linked to stability patterns, because the asynchrony of local population dynamics leads to more stable dynamics (low biomass Table 1: Approximative relative increase in predation risk between low-risk and high-risk environments (equivalent to the asymmetry of interaction strength γ in Figure 1). See [START_REF] Gorini | Habitat heterogeneity and mammalian predator-prey interactions: Predator-prey interactions in a spatial world[END_REF] CV) at higher scales due to compensatory dynamics [START_REF] Loreau | Biodiversity as spatial insurance in heterogeneous landscapes[END_REF][START_REF] Gonzalez | The causes and consequences of compensatory dynamics in ecological communities[END_REF][START_REF] Loreau | Biodiversity and ecosystem stability: A synthesis of underlying mechanisms[END_REF][START_REF] Wilcox | Asynchrony among local communities stabilises ecosystem function of metacommunities[END_REF], subsequent studies suggested that increased asymmetry does not necessarily leads to increased stability. For example [START_REF] Ruokolainen | The roles of spatial heterogeneity and adaptive movement in stabilizing (or destabilizing) simple metacommunities[END_REF] presented a model in which biomass fluctuations can become more variable with increasing asynchrony. Hence, the relationship between asymmetry and stability is not trivial and the mechanisms governing asynchrony through the difference in energy flow between the fast ans slow channels are not well understood. To fill this gap, we propose to consider the effects of asymmetry from the metacommunity perspective since recent theoretical studies were able to accurately explain the synchrony and stability patterns in metacommunities (Quévreux et al., 2021a,b). Metacommunities embody the spatial dimension of interaction networks: they consist of distant patches connected by the dispersal of the organisms living in each patch [START_REF] Leibold | The metacommunity concept: A framework for multi-scale community ecology[END_REF]Leibold and Chase, 2017). The metacommunity framework is particularly suitable to represent the spatial heterogeneity observed in ecosystems because each community has its own characteristics such as biomass turnover. Following Rooney et al.'s (2006) model, many studies implemented spatial heterogeneity through the asymmetry of interaction strength and/or resource supply to manipulate the difference in biomass turnover between the energy channels hosted by each patch [START_REF] Goldwyn | Small heterogeneity has large effects on synchronization of ecological oscillators[END_REF][START_REF] Ruokolainen | The roles of spatial heterogeneity and adaptive movement in stabilizing (or destabilizing) simple metacommunities[END_REF][START_REF] Anderson | Body size dependent dispersal influences stability in heterogeneous metacommunities[END_REF]. In particular, interaction strength is key in community dynamics because it governs food web structure, stability [START_REF] Neutel | Stability in real food webs: Weak links in long loops[END_REF] and biomass distribution (Barbier and Loreau, 2019) by simultaneously determining predator growth and prey mortality. Therefore, its significant variations observed in nature, often reported as predation risk by prey in field studies (Table 1), should lead to dramatic variations in community functioning across space.

In addition to the asymmetry of interaction strength, Rooney et al. (2006) highlighted the im-portance of mobile predators coupling two different energy channels, a keystone role in ecosystem functioning largely reported by empirical studies [START_REF] Schindler | Habitat coupling in lake ecosystems[END_REF][START_REF] Vadeboncoeur | Effects of multi-chain omnivory on the strength of trophic control in lakes[END_REF][START_REF] Schmitz | Perturbation and abrupt shift in trophic control of biodiversity and productivity: Perturbation and regime shift[END_REF][START_REF] Olff | Parallel ecological networks in ecosystems[END_REF][START_REF] Dolson | Lake morphometry predicts the degree of habitat coupling by a mobile predator[END_REF][START_REF] Schmitz | Predator control of ecosystem nutrient dynamics[END_REF]. In Rooney et al.'s (2006) model, the perturbation of the mobile predator leads to an asynchronous response of prey, which stabilises the food web. However, Quévreux et al. (2021a) showed that the perturbation and dispersal of particular trophic levels govern synchrony and stability in symmetric metacommunities. In asymmetric metacommunities, the perturbation of particular patches should lead to different synchrony and stability patterns at the metacommunity scale because of the different dynamics in each patch in response to perturbations. In parallel to the keystone role of mobile predators, keystone communities (sensus [START_REF] Mouquet | Extending the concept of keystone species to communities and ecosystems (B. Blasius[END_REF], which are equivalent to keystone patches), should have a major influence on synchrony and stability patterns. Therefore, we expect that asymmetry is not a generic stabilising factor, as claimed by Rooney et al. (2006), but strongly depends on which patch is perturbed according to its characteristics. To explore this statement, we consider a simple metacommunity model of two patches hosting a predator-prey couple and with asymmetric interaction strength and resource supply. The stability of the metacommunity is assessed by the response at different scales (e.g. CV of the biomass of a species at the local and regional scales) when prey receive stochastic perturbations in one of the two patches. We show contrasting effects of asymmetry: perturbing prey in the fast patch (equivalent to the fast channel defined by [START_REF] Rooney | Integrating food web diversity, structure and stability[END_REF] promotes prey synchrony and decreases predator stability at the metapopulation scale while perturbing the slow patch has the opposite effects.

Positive effect of the prey on its predator

Negative effect of the predator on its prey Self-regulation Growth rate patch #1 fast patch patch #2 slow patch asymmetry of interaction strength asymmetry of resource supply stochastic perturbations dispersal Figure 1: The metacommunity model consists of two patches, each sustaining a predator-prey couple linked by predators, which disperse at a very high scaled rate d 2 . Prey grow at a rate g/D and have a positive effect εa on predators, while predators have a negative effect ma on prey. Each species population is also limited by self-regulation D (negative intraspecific interactions). Spatial heterogeneity is embodied by the asymmetry of resource supply (green) and the interaction strength (purple), which are higher in patch #1 by factors ω and γ, respectively. Consistent with Rooney et al. (2006), patch #1 is called the fast patch, and patch #2 is called the slow patch. Prey receive stochastic perturbations either in patch #1 or in patch #2 (red arrows).

1 D dB (1) 1 dt = B (1) 1 ω g D -B (1) 1 -γmaB (1) 2 (1a) 1 mD dB (1) 2 dt = B (1) 2 -B (1)
2 + γϵaB

(1) 1 intra-patch dynamics + d 2 B (2) 2 -B (1) 2 dispersal (1b)

Methods

Metacommunity model

We use the model proposed by Quévreux et al. (2021a) based on the food chain model developed by Barbier and Loreau (2019). The model consists of two patches that each sustain a food chain with Lotka-Volterra predator-prey interactions (equations (1a) and (1b)).

B (k) 1 and B (k) 2
are the biomasses of prey and predators, respectively, in patch k. Prey have a positive effect εa on predators (ε is the conversion efficiency and a is the interspecific interaction rate relative to intraspecific interactions), and predators have a negative effect ma on prey (m is the predator to prey metabolic rate ratio)(Figure 1). The time scale of the system is rescaled by the metabolic rate of prey, and biologic rates of each species i are rescaled by its intraspecific interaction rate D i . Therefore, we obtain the relative growth rate g/D and the scaled dispersal rate d i . Considering scaled parameters and aggregated parameters (εa and ma) enables us to explore a wide range of ecological situations. We refer to Appendix S1-1 for a detailed description of the food chain model and analysis methods. All the necessary information to fully understand our results are in the main text and the supplementary information only serves to give additional technical elements to fully reproduce our work and proofs of the robustness of our results. Parameters and their values are summarised in Table 2.

We reproduce the two main features of Rooney et al.'s (2006) model. First, predators disperse at a very high scaled rate d 2 = 10 6 , while prey are immobile (d 1 = 0), and strongly couple the two patches. Slightly mobile prey (0 < d 1 ≪ d 2 ) should not change the results because Quévreux et al. (2021a) showed that the species for which dispersal has the strongest influence drives the coupling between the two patches. Second, resource supply and interaction strength are asymmetric between patches since they are higher in patch #1 by factors γ and ω respectively (Figure 1). Patch #1 corresponds to the fast energy channel, in which biomass has a high turnover, while patch #2 corresponds to the slow channel. Therefore, we call patch #1 the fast patch and patch #2 the slow patch. We set γ = ω to ensure species persistence over the entire range of parameters (see Figure S2-6 in the supporting information) but varying them independently does not qualitatively change the results (see Figure S2-14 in the supporting information). In the following, we only refer to γ for the sake of simplicity and only consider γ ≥ 1 because γ ≤ 1 just swaps the roles of patches #1 and #2.

Table 2: Table of parameters. σ i is set very small to keep the system in the vicinity of equilibrium. More combinations of εa and ma are tested in the supporting information. d 2 is set very high to emphasise the high mobility of predators and their ability to couple prey populations. ω is set equal to γ. r = 0 removes the energetic limitations of the food chain and makes interactions the dominant factors determining biomass distribution and stability patterns, as in Barbier and Loreau (2019)

Response to stochastic perturbations

We use the same methods as Quévreux et al. (2021a) to study the response of metacommunities to stochastic perturbations. Indeed, recent studies advocate for the use of the temporal variability of biomass (Haegeman et al., 2016;Arnoldi et al., 2018), which is measured by the coefficient of variation (CV), and can be easily measured experimentally. In addition, Wang andLoreau (2014, 2016), Wang et al. (2019), andJarillo et al. (2022) showed that CVs scale up from local populations to community, regional and metacommunity levels, therefore providing a comparison of stability at different scales. Here, we provide only a brief description of the main concepts, but a thorough description is available in Appendix S1.

Prey in the fast or slow channel receive stochastic perturbations that are represented by equation (2).

dB i = f i (B 1 , ..., B S )dt Deterministic + σ i B * i dW i Perturbation (2) f i (B 1 , ..., B S )
represents the deterministic part of the dynamics of species i, as described by equations (1a) and (1b)). Stochastic perturbations are defined by their standard deviation σ i and dW i , a white noise term with a mean of 0 and variance of 1. Perturbations also scale with the square root of the biomass at equilibrium B * i of the perturbed population. Such scaling makes the perturbations similar to demographic stochasticity (from birth-death processes) that evenly affect each species regardless of abundance (Arnoldi et al., 2019). In other words, the ratio of mean species biomass variance to perturbation variance is roughly independent of biomass, which disentangles the effect of asymmetry on perturbation transmission from its effect on species abundance. Therefore, for different per-turbations affecting different species with the same value of standard deviation σ i , we generate a similar variance at the metacommunity scale regardless the abundance of the perturbed species and excite the entire metacommunity with the same intensity (see Figure S2-3 in the supporting information).

In the following, we assess the temporal variability of the biomass of each population induced by stochastic perturbations affecting the metacommunity. Therefore, we linearise the system in the vicinity of equilibrium to obtain equation (3) where

X i = B i -B * i is the deviation from equilibrium. d - → X dt = J - → X + T - → E (3)
J is the Jacobian matrix, which represents the linearised direct effects between populations in the vicinity of equilibrium, and T defines how the perturbations E i = σ i dW i apply to the system (i.e., which species they affect and how they scale with biomass, where T is a diagonal matrix whose terms are T ii = B * i ). Because the system is at steady state, the stationary variance-covariance matrix C * of species biomasses (variance-covariance matrix of -→ X , see the demonstration in Appendix S1-5) can be calculated from the variance-covariance matrix of perturbations V E (variance-covariance matrix of -→ E ) by solving the Lyapunov equation (4) [START_REF] Arnold | Stochastic differential equations: Theory and applications[END_REF]Wang et al., 2015;Arnoldi et al., 2016;Quévreux et al., 2021a).

JC * + C * J ⊤ + T V E T ⊤ = 0 (4)
The expressions for V E and T and the method to solve the Lyapunov equation are detailed in Appendix S1-6. From the variance-covariance matrix C * , we compute the coefficient of correlation of the biomass dynamics between the two populations of each species (see equation ( 22) in the supporting information) and we measure the stability with the coefficient of variation (CV) of the biomass. In addition, biomass CVs can be measured at different scales: population scale (e.g., biomass CV of prey in patch #1), metapopulation scale (e.g., CV of the total biomass of prey) and metacommunity scale (e.g., CV of the total biomass of predator and prey put together) to assess the effects of asymmetry at local and regional scales (Figure 3A and see Appendix S1-7). Finally, we quantify the synchrony of the dynamics of the different populations with the coefficient of correlation, which is also computed from the variance-covariance matrix C * (Appendix S1-7).

Results

Effects on stability

We describe how the asymmetry of interaction strength γ shapes metacommunity stability at different scales. Since predators have a very high scaled dispersal rate (d 2 = 10 6 ), their populations are perfectly correlated and display the same dynamics. Our main result is that prey become more correlated when they are perturbed in patch #1 (fast channel in which γ > 1), while they become more anticorrelated when they are perturbed in patch #2 (Figure 2). Increasing γ amplifies the difference in correlation between these two scenarios, and this pattern qualitatively holds for various combinations of the physiological and ecological parameters εa and ma (see Figure S2-9 in the supporting information).

Increasing the asymmetry of interaction strength γ has contrasting effects on biomass CV at different scales as well (Figure 3A). At the population scale, it increases the biomass CV of each population when prey are perturbed in the fast patch (Figure 3B). When prey are perturbed in the slow patch, increasing γ slightly alters the biomass CV of prey in patch #1, increases the biomass CV of prey in patch #2 and decreases the biomass CV of predators. This discrepancy can be attributed to the strong effect of γ on prey biomass in patch #2 (Figure 4): prey biomass strongly decreases with γ in patch #2, which increases their biomass CV.

At the metapopulation scale, the asymmetry of interaction strength γ increases the biomass CV of prey in both scenarios of perturbation (Figure 3C). However, this result is not true for all values of εa and ma (Figure S2-10A in the supporting information) because of the various responses of prey biomass to γ among patches (Figure S2-8A in the supporting information). The biomass CV of predators is higher when prey are perturbed in the fast patch (patch #1) compared to the case in which prey are perturbed in the slow patch (#2) (Figure 3C), which is consistent for all values of εa and Figure 2: Spatial correlation between the populations of each species depending on asymmetry of interaction strength γ when predators disperse and prey are perturbed in patch #1 or #2. The dashed line emphasises the value of the correlation of prey populations without asymmetry (γ = 1). Note that the curves for predators overlap because their high dispersal that perfectly correlates their dynamics regardless of the perturbed patch. ma (see Figure S2-10A in the supporting information).

Finally, stability at the metacommunity scale depends on the distribution of biomass and CV among species. In our particular case (εa = 1 and ma = 1), predators have the largest total biomass (Figure 4) and drive the biomass CV at the metacommunity scale for low values of asymmetry of interaction strength γ (Figure 3D). For high values of γ, when prey are perturbed in patch #2, the CV of total biomass increases with γ because it is driven by prey in patch #2, whose biomass CV is much higher than the biomass CV of predators, which compensates for their lower biomass. Other values of εa and ma lead to other distributions of biomass and CV among species, which can make prey to drive the stability at the metacommunity scale (see Fig-

Underlying mechanisms

To unveil the mechanisms governing the stability of heterogeneous metacommunities, we look deeper into the dynamics after a pulse perturbation (Figure 5A) and explain them with the direct effects between species quantified by the Jacobian matrix (see equation ( 3)). When the perturbation of prey occurs in patch #1, the strong direct effect of prey on predators (and vice versa) in patch #1 due to γ (Figure 5B) leads to a strong response of predators both patches and correlate their dynamics, which explains why asymmetry of interaction strength γ increases prey correlation when prey are perturbed in patch #1.

When the perturbation of prey occurs in patch #2, the weak direct effect of prey on predators (Figure 5B) leads to a small response of predators 3

⃝.

In turn, the very low direct effect of predators on prey in patch #2 does not allow perturbations to ripple back to patch #2 where prey slowly respond from the initial perturbation 4 ⃝ (Figure 5B). This slow response is emphasised by the source-sink dynamics in the metacommunity (Figure 4 and Figure S2-5 in the supporting information), which leads to a lower biomass of prey in patch #2 compared to a metacommunity without dispersal, therefore decreasing biomass flows in patch #2 and its response speed. This difference in response speed between patches #1 and #2 leads to the anti-correlation of prey populations because it increases the time interval in which they have opposite variations: an increase in the biomass of prey in patch #1, which follows the initial decrease due to predation, and a slow decrease in prey biomass in patch #2.

Discussion

We have shown that the asymmetry of interaction strength has contrasting effects on stability depending on which patch is perturbed. Perturbing prey in the fast patch (in which interaction strength is the highest) tends to synchronise the dynamics of prey biomass and increases the temporal variability of predator dynamics at the metapopulation scale, while perturbing prey in the slow patch decreases the synchrony of prey dynamics and the temporal variability predator dynamics. This discrepancy between the responses is due to asymmetric transmission of perturbations within each patch, itself caused by different biomass distributions and interaction strengths between the two patches. Perturbations are strongly transmitted from the fast patch to the slow patch, while the reverse transmission is weak. Consequently, the fast patch drives the dynamics of the metacommunity and synchronises prey dynamics, while the slow patch does not, and the almost independent respond from perturbation in each patch decreases the synchrony of prey dynamics.

Stability in a heterogeneous world

Our results show that spatial heterogeneity, which is represented by the asymmetry of interaction strength and resource supply as in Rooney et al. (2006), generates mechanisms that alter local and regional dynamics, which deeply changes the synchrony of population dynamics and the stability of the metacommunity at different scales. Quévreux et al. (2021a) showed that, in a homogenous metacommunity, the spatial correlations between patches can be obtained from the withinpatch correlations, the dispersing species making the link between the two (see Figure S2-28 for a summary of the results of Quévreux et al., 2021a). In other words, knowledge of the dynamics at the local scale is enough to understand the stability pattern at the metacommunity scale. In a heterogeneous metacommunity, a similar approach does not work because patches do not contribute equally to the dynamics. In a heterogeneous metacommunity, a similar approach does not work because patches do not contribute equally to the dynamics. In particular, a patch with fast energy flow can have an overwhelming impact (see Figure S2-7 in the supporting information). Clearly, the dynamics at the metacommunity scale cannot be assessed by the dynamics at the local scale, as in Quévreux et al. (2021a), and they are an emergent property resulting from the tight interplay between the strength of perturbation transmission in each patch. Rooney et al. (2006) verbally explained that each patch has dynamics with different speeds: the fast channel (with higher interaction rates and resource supplies) enables a quick response after a perturbation while the slow channel dampens the dynamics in the long term and prevents the system from overshooting. By considering the stability at different scales, our results contrast with this explanation. On the one hand, the asynchrony of prey dynamics, when they are perturbed in patch #2 (Figure 2), stabilises the dynamics of predators because their resource supplies are asynchronous. On the other hand, the dynamics of prey at the metapopulation scale are not stabilised by their asynchrony (Figure 3C) because of the low local stability in patch #2 (Figure 3B), which decreases the overall stability of prey. The potential stabilising effect of asymmetry depends both on the perturbed patch and the considered trophic level. Therefore, the overall stability at the metacommunity scale is governed by the relative contributions of the various populations in response to local perturbations, and asymmetry per se does not have a stabilising effect.

Our description of the mechanisms underlying the apparent stabilising effects of spatial heterogeneity should enlighten the results of previous theoretical studies. [START_REF] Goldwyn | Small heterogeneity has large effects on synchronization of ecological oscillators[END_REF] and [START_REF] Ruokolainen | The roles of spatial heterogeneity and adaptive movement in stabilizing (or destabilizing) simple metacommunities[END_REF] found that the asymmetry of interaction rate leads to asynchrony by generating out-of-phase dynamics in a system with endogenous oscillations. In particular, [START_REF] Ruokolainen | The roles of spatial heterogeneity and adaptive movement in stabilizing (or destabilizing) simple metacommunities[END_REF] found a U-shaped relationship: for moderate asymmetry, the spatial asynchrony of predator and prey populations is maximal, which leads to optimum stability at the metacommunity scale. Our results suggest that moderate asymmetry would alter the phase of the oscillations in each patch while keeping the amplitude of oscillations equivalent, therefore promoting asynchrony. Conversely, a strong asymmetry would increase the imbalance between oscillation amplitude and enable the fast patch to take over the slow patch, which would bring back synchrony. However, their results rely on phase-locking [START_REF] Jansen | Phase locking: Another cause of synchronicity in predator-prey systems[END_REF][START_REF] Lloyd | Synchronicity, chaos and population cycles: Spatial coherence in an uncertain world[END_REF][START_REF] Goldwyn | When can dispersal synchronize populations?[END_REF][START_REF] Vasseur | Phase-locking and environmental fluctuations generate synchrony in a predator-prey community[END_REF], which is the coupling of the phase of oscillators embodied by predator-prey pairs in each patch. Although our results provide interesting insight into metacommunity dynamics, they can-not grasp the fine mechanisms underlying nonlinear phenomena such as phase-locking and further studies are needed to identify these mechanisms.

Generality of the effects of asymmetry on stability

Our main results is that asymmetry is stabilising when the slow patch is perturbed, while it is destabilising when the fast patch is perturbed. This result is strikingly robust to several deviations from the original model we have described. First, we show that the described mechanisms are valid for a wide range of ecological and physiological parameters leading to various distributions of biomass among predators and prey (see Figures S2-8 andS2-9 in the supporting information). In addition, we observe the same results for longer food chains as long as prey populations are directly coupled by the dispersing predator (see Figures S2-18 and S2-19 in the supporting information). Currently, we do not identify a clear pattern for species lower in the food chain over a wide range of ecological and physiological parameters but further studies are needed to investigate the potential indirect effects propagating across the food chain. Second, the mechanisms are not restricted to prey populations coupled by a mobile predator but also apply to predator populations coupled by a mobile prey (see Figure S2-22 in the supporting information). Therefore, we anticipate that mobile predators are not the only major drivers of synchrony and stability in ecosystems [START_REF] Schindler | Habitat coupling in lake ecosystems[END_REF][START_REF] Vadeboncoeur | Effects of multi-chain omnivory on the strength of trophic control in lakes[END_REF][START_REF] Dolson | Lake morphometry predicts the degree of habitat coupling by a mobile predator[END_REF][START_REF] Olff | Parallel ecological networks in ecosystems[END_REF][START_REF] Rooney | Integrating food web diversity, structure and stability[END_REF], and resource species may also have an equivalent impact. Taken together, these two points strongly suggest that the mechanisms underlying stability and synchrony in response to perturbations should be general to metacommunities regardless of the ecological parameters, biomass distribution and dispersal among species.

Spatial heterogeneity has often been presented as a generic condition generating mechanisms ensuring stability, but our results contradict this statement. The models focusing on the asymmetric feeding of consumers on different energy channels or different patches showed that it promotes the existence of stable equilibria [START_REF] Mccann | Weak trophic interactions and the balance of nature[END_REF], greater asymptotic resilience (Rooney et al., 2006), asynchrony of prey in response to predator perturbation (Rooney et al., 2006) and out-of-phase limit cycles [START_REF] Goldwyn | Small heterogeneity has large effects on synchronization of ecological oscillators[END_REF][START_REF] Ruokolainen | The roles of spatial heterogeneity and adaptive movement in stabilizing (or destabilizing) simple metacommunities[END_REF]. All these studies considered measures of stability aiming to capture the general stability properties of metacommunities and miss the targeted effects of perturbations as we explained earlier. Although asymmetry does not necessarily promote stability, our results show that general mechanisms drive the response of metacommunities to localised perturbations, therefore providing a valuable framework to assess the response of ecosystems to localised perturbations due to human activity. Additionally, these mechanisms enable us to understand the effect of environmental perturbations affecting all patches. As demonstrated by Arnoldi et al. (2019), environmental perturbations affect abundant populations the most, which is the prey population in the fast patch in our case (see Figure S2-25 in the supporting information). Therefore, we anticipate that the fast patch will govern the dynamics of metacommunities in which all populations are perturbed (see Figure S2-27 in the supporting information).

Implications for conservation

The metacommunity framework has long been used in conservation ecology [START_REF] Johnson | Urbanization and wetland communities: Applying metacommunity theory to understand the local and landscape effects (M. Cadotte[END_REF][START_REF] Schiesari | Towards an applied metaecology[END_REF][START_REF] Patrick | The application of metacommunity theory to the management of riverine ecosystems[END_REF]. Conservation efforts are usually concentrated on particular locations and useful management must consider the ecological processes acting at the landscape scale [START_REF] Van Teeffelen | Species in a dynamic world: Consequences of habitat network dynamics on conservation planning[END_REF][START_REF] Chase | Biodiversity conservation through the lens of metacommunity ecology[END_REF]. For instance, spatial heterogeneity is key to ensuring species coexistence and diversity at the regional scale, which ultimately provides important ecosystem services in agricultural landscapes [START_REF] Bennett | Properties of land mosaics: Implications for nature conservation in agricultural environments[END_REF]. A large corpus of theoretical studies explored the local response of communities in a landscape receiving perturbations [START_REF] Mouquet | Extinction debt in source-sink metacommunities (T. Romanuk[END_REF][START_REF] Economo | Biodiversity conservation in metacommunity networks: Linking pattern and persistence[END_REF][START_REF] Holyoak | Integrating disturbance, seasonality, multi-year temporal dynamics, and dormancy into the dynamics and conservation of metacommunities[END_REF][START_REF] Jacquet | Temporal variation of patch connectivity determines biodiversity recovery from recurrent disturbances[END_REF]. However, these studies focus on extinction events recovered by dispersal events in a patch dynamics framework, and little is known about the effect of moderate or small perturbations. In this context, the present study provides valuable insight into fine-scale dynamics in response to perturbations.

Our results show that species interactions are a major driver of synchrony in heterogeneous metacommunities. Even if the species of interest does not disperse significantly, the synchrony of the dynamics of its different populations can strongly depend on the interactions with another species with a higher dispersal across the landscape. For instance, [START_REF] Howeth | Predation inhibits the positive effect of dispersal on intraspecific and interspecific synchrony in pond metacommunities[END_REF] showed that predatory fish promote the asynchrony of oscillating populations of zooplankton in a mesocosm experiment. Therefore, species endorsing this role are called "mobile link organisms" [START_REF] Lundberg | Mobile link organisms and ecosystem functioning: Implications for ecosystem resilience and management[END_REF] and are particularly targeted by conservation policy because they have major impacts on community dynamics and ecosystem functioning [START_REF] Soulé | Strongly interacting species: Conservation policy, management, and ethics[END_REF][START_REF] Brodie | Ecological function analysis: Incorporating species roles into conservation[END_REF]. Such a species can be considered keystone species [START_REF] Mills | The keystone-species con-cept in ecology and conservation[END_REF] and must be clearly identified to properly manage the conservation of the other interaction species. However, our results show that mobile link organisms are not the only driver of metacommunity stability, and the patch being perturbed also has a major impact. The concept of a keystone community, defined by [START_REF] Mouquet | Extending the concept of keystone species to communities and ecosystems (B. Blasius[END_REF] for communities whose destruction causes species extinction or a decrease in biomass production, can be applied to better assess the stability of metacommunities. Keystone communities are usually identified as those patches that are strongly connected to other patches in the spatial network [START_REF] Resetarits | Testing the keystone community concept: Effects of landscape, patch removal, and environment on metacommunity structure[END_REF], but our results suggest that the dynamical properties of each patch can be important as well. For instance, the fast patch can be identified as a keystone patch because of its ability to synchronise the dynamics of the other patches. Therefore, identifying the communities living in fast and slow patches should be key for conservation management aiming to mitigate the effects of perturbations.

According to our results, mitigating the effects of perturbations affecting the patch in which interaction strength is the highest is critical to avoid the synchrony of prey dynamics (Figure 2) and ensure predator stability (Figure 3C). Then, the patch in which the interaction strength between the species of interest and the mobile link organism is the highest must be identified. Conservation policies usually target preserved areas because they are characterised by high species richness but identifying them as fast or slow patches is not trivial. Urban ecology is a relevant example because many species dwell in cities and less anthropised ecosystems (e.g., agricultural and natural landscapes). Urban areas can be considered fast patches because of the abundance of resources (parameters ω in our model) for opportunistic species, but they can also be considered slow patches because of the reduced predation pressure (parameter γ in our model), cities acting as safe spaces (see [START_REF] Shochat | From patterns to emerging processes in mechanistic urban ecology[END_REF] and [START_REF] Shochat | Invasion, competition, and biodiversity loss in urban ecosystems[END_REF] for review). Typically, birds and rodents can find plenty of food due to human wastes, public parks and feeding while experiencing less predation [START_REF] Rebolo-Ifrán | Urban conservation hotspots: Predation release allows the grassland-specialist burrowing owl to perform better in the city[END_REF]. Therefore, focusing conservation efforts on urban areas to mitigate the perturbations affecting their ecosystem may be as important as protecting wild areas to protect species at the metapopulation scale.

Conclusion

Asymmetry of interaction strength, and spatial heterogeneity in general, is not stabilising factor per se because perturbing prey in the fast patch leads to the synchrony of the dynamics of prey populations and increases the temporal variability of the mobile predator linking the two patches. Therefore, the response of metacommunities to perturbations is strongly context dependent, i.e., a good knowledge of the characteristics of each patch relative to each other is required to assess stability at the metacommunity scale. Based on our findings, we advocate for conservation efforts to target key patches not only according to species richness or biomass density but also according to the distribution of interaction strength across the metacommunity.

S1 Complementary material and methods

S1-1 Model description

The model has been originally developed by Barbier and Loreau (2019), who considered a food chain model with a simple metabolic parametrisation. Their model corresponds to the "intra-patch dynamics" part of equations ( 5a) and (5b) to which we graft a dispersal term to consider a metacommunity with two patches.

dB (1) 1 dt = B (1) 1 (ω 1 g 1 -D 1 B (1) 1 -γ 1 α 2,1 B (1) 2 ) + δ 1 (B (2) 1 -B (1) 1 ) (5a) dB (1) i dt = B (1) i (-r i -D i B (1) i + γ 1 ϵα i,i-1 B (1) i-1 -γ 1 α i+1,i B (1) i+1 ) intra-patch dynamics + δ i (B (2) i -B (1) i ) dispersal (5b) B (1) i
is the biomass of trophic level i in the patch #1, ϵ is the biomass conversion efficiency and α i,j is the interaction strength between consumer i and prey j. Species i disperses between the two patches at rate δ i . The density independent net growth rate of primary producers g i in equations (5a), the mortality rate of consumers r i in equations (5b) and the density dependent mortality rate D i scale with species metabolic rates m i as biological rates are linked to energy expenditure.

g 1 = m 1 g r i = m i r D i = m i D (6)
In order to get a broad range of possible responses, we assume the predator-prey metabolic rate ratio m and the interaction strength to self-regulation ratio a to be constant. These ratios capture the relations between parameters and trophic levels. This enables us to consider contrasting situations while keeping the model as simple as possible.

m = m i+1 m i a = α i,i-1 D i d i = δ i D i (7)
Varying m leads to food chains where predators have faster or slower biomass dynamics than their prey and varying a leads to food chains where interspecific interactions prevail or not compared with intraspecific interactions. As all biological rates are rescaled by D i , we also define d i , the dispersal rate relative to self-regulation (referred as scaled dispersal rate in the rest of the study), in order to keep the values of the dispersal rate relative to the other biological rates consistent across trophic levels. Finally, the time scale of the system is defined by setting the metabolic rate of the primary producer m 1 to unity. Thus, we can transform equations (5a) and (5b) into:

1 D dB (1) 1 dt = B (1) 1 (ω g D -B (1) 1 -γmaB (1) 2 ) + d 1 (B (2) 1 -B (1) 1 ) (8a) 1 m i-1 D dB (1) i dt = B (1) i (- r D -B (1) i + γϵaB (1) i-1 -γmaB (1) i+1 )
intra-patch dynamics

+ d i (B (2) i -B (1) i ) dispersal (8b)
Thus, ϵa and ma defines the positive effect of the prey on its predator and the negative effect of the predator on its prey, respectively. These two synthetic parameters define the overall behaviour of the food chain and will be varied over the interval [0.1, 10] to consider a broad range of possible responses. Finally, the mortality rate is set to zero (r = 0) to remove the energetic limitations of the food chain and make interactions the dominant factors determining biomass distribution and stability patterns, as in Barbier and Loreau (2019).

S1-2 Biomass at equilibrium when top predators populations are perfectly coupled

The system can be easily solved if we consider the total population of top predator instead of two populations connected by dispersal. Since the two populations are perfectly coupled by dispersal, top predator i biomass is constant across patches and we have B (2) * i = 0.5B tot * i . Thus we have the following system at equilibrium for the two top species (the equations for the other species are the same as the symmetric case):

0 = - r D -B * tot i + εa γB * (1) i-1 + B * (2) i-1 (9a) 0 = - r D -B * (1) i-1 + γεaB * (1) i-2 -γma B * tot i 2 (9b) 0 = - r D -B * (2) i-1 + εaB * (2) i-2 -ma B * tot i 2 (9c)
Which translates into the following matrix equation:

                -1 -γma γεa -1 -γma . . . . . . . . . (0) γεa -1 0 -γma/2 0 -1 -ma εa -1 -ma (0) . . . . . . . . . εa -1 -ma/2 γεa εa -1       ωg/D -r/D . . . -r/D g/D -r/D . . . -r/D -r/D                 = 0 (10)

S1-3 Biomass at equilibrium when basal species populations are perfectly coupled

In the same way, we have:

0 = g D -B * tot 1 -ma γB * (1) 2 + B * (2) 2 (11a) 0 = - r D -B * (1) 2 + γεa B * tot 1 2 -γmaB * (1) 3 (11b) 0 = - r D -B * (2) 2 + εa B * tot 1 2 -maB * (2) 3 (11c)
Which translates into the following matrix equation:

            -1 -γma -ma γεa/2 -1 -γma (0) . . . . . . . . . γεa -1 0 εa/2 0 -1 -ma (0) . . . . . . . . . εa -1                          B * tot 1 B * (1) 2 . . . B * (1) i B * (2) 2 . . . B * (2) i              +             (ω + 1)g/D -r/D . . . -r/D -r/D . . . -r/D             = 0 (12)

S1-4 Linearisation of the system

The system of equations ( 1a) and (1b) can be linearised in the vicinity of equilibrium:

dB i dt = f i (B * 1 , ..., B * S ) = 0 + S j=1 ∂f i ∂B j B * (B j -B * j ) (13)
Thus, by setting X i = B i -B * i the deviation from equilibrium, we have:

dX i dt = S j=1 J ij X j (14)
Then, we can consider small perturbations defined by -→ E whose effects on -→ X are defined by the matrix T (Arnoldi et al., 2016). We get the linearised version of equation ( 2

): d - → X dt = J - → X + T - → E (15)
The elements of -→ E are defined by stochastic perturbations E i = σ i dW i with σ i their standard deviation and dW i a white noise term with mean 0 and variance 1. In our model, each species i in each patch k can receive demographic perturbations scaling with the square root of their biomass at equilibrium. Thus, -→ E contains the white noise term σ

(k) i dW (k) i
for each population of each species, T is a diagonal matrix whose terms are B * (k) i and the matrix product T -→ E results in the product of the white noise and the biomass scaling as in equation ( 2) in the main text.

S1-5 Demonstration of the Lyapunov equation

The following demonstration of the Lypunov equation has been taken from [START_REF] Oku | On the covariance matrix of the stationary distribution of a noisy dynamical system[END_REF]. The continuous-time dynamics from equation ( 15) can be converted to a discrete-time dynamics by using Euler-Maruyama method:

- → X t+∆t = - → X t + ∆tJ - → X t + √ ∆tT - → E t (16) C * = E[ - → X t - → X ⊤ t ] (the expected value of the product - → X t and its transpose - → X ⊤ t )
is the stationary variancecovariance matrix of the system, therefore dC * /dt = 0. We also have the following relation:

dC * dt = lim ∆t→0 E[ - → X t+∆t - → X ⊤ t+∆t ] -E[ - → X t - → X ⊤ t ] ∆t = lim ∆t→0 E[( - → X t + ∆tJ - → X t + √ ∆tT - → E t )( - → X t + ∆tJ - → X t + √ ∆tT - → E t ) ⊤ ] -E[ - → X t - → X ⊤ t ] ∆t = lim ∆t→0 ∆tE[ - → X t - → X ⊤ t J ⊤ ] + ∆tE[J - → X t - → X ⊤ t ] + ∆t 2 E[J - → X t - → X ⊤ t J ⊤ ] + ∆tE(T - → E t - → E ⊤ t T ⊤ ) ∆t = E[ - → X t - → X ⊤ t J ⊤ ] + E[J - → X t - → X ⊤ t ] + E[T - → E t - → E ⊤ t T ⊤ ] = C * J ⊤ + JC * + T V E T ⊤ = 0 (17) Because E[ - → X t ] = 0, E[ - → E t ] = 0, E[ - → X t - → E ⊤ t ] = 0, E[ - → E t - → X ⊤ t ] = 0 and V E = E[ - → E t - → E ⊤ t ]
the variance-covariance matrix of stochastic perturbations.

S1-6 Resolution of the Lyapunov equation

In the vicinity of equilibrium, the Lyapunov equation links the variance-covariance matrix V E of the perturbation vector -→ E to the variance-covariance matrix C * of species biomasses (see the appendix of Wang et al. (2015) for more details on the Lyaponov equation).

JC * + C * J ⊤ + T V E T ⊤ = 0 (18)
The diagonal elements of V E are equal to σ 2 i (variance of the white noises) and the non-diagonal elements are equal to zero because perturbations are independent. ⊤ is the transpose operator. C * can be calculated using a Kronecker product [START_REF] Nip | Direct numerical solution of algebraic Lyapunov equations for largescale systems using Quantized Tensor Trains[END_REF]. The Kronecker product of an m × n matrix A and a p × q matrix B denoted A ⊗ B is the mp × nq block matrix given by:

A ⊗ B =    a 11 B • • • a 1n B . . . . . . . . . a m1 B • • • a mn B   
We define C * s and (T V E T ⊤ ) s the vectors stacking the columns of C * and T V E T ⊤ respectively. Thus, equation ( 18) can be rewrite as:

(J ⊗ I + I ⊗ J)C * s = -(T V E T ⊤ ) s C * s = -(J ⊗ I + I ⊗ J) -1 (T V E T ⊤ ) s (19)

S1-7 Coefficient of variation and correlation

Our different metrics of stability can be easily computed from the elements of the variance-covariance matrix C * defined by elements w i (k) j (ℓ) that are the covariance between species i in patch k and species j in patch ℓ.

C * =              w 1 (1) 1 (1) • • • w 1 (1) S (1) • • • w 1 (n) 1 (n) • • • w 1 (n) S (n) . . . . . . . . . • • • . . . . . . . . . w S (1) 1 (1) • • • w S (1) S (1) • • • w S (n) 1 (n) • • • w S (n) S (n) . . . . . . . . . . . . . . . . . . . . . w 1 (n) 1 (1) • • • w 1 (n) S (1) • • • w 1 (n) 1 (n) • • • w 1 (n) S (n) . . . . . . . . . • • • . . . . . . . . . w S (n) 1 (1) • • • w S (n) S (1) • • • w S (n) 1 (n) • • • w S (n) S (n)              (20)
The temporal variability of the metacommunity is assessed with the coefficient of variation (CV) of biomass at different scales: population scale CV

(k) i
, which is the biomass CV of species i in patch k, metapopulation scale CV i , which is the biomass CV of the total biomass of species i across patches and metacommunity scale CV M C , which is the total biomass of the entire metacommunity (Wang and Loreau, 2014;Wang et al., 2019;Jarillo et al., 2022).

CV (k) i = √ w i (k) i (k) µ (k) i CV i = kℓ w i (k) j (ℓ) k B * (k) i CV M C = ijkℓ w i (k) j (ℓ) ik B * (k) i (21)
Figure S1-1: Elements of the variance-covariance matrix C * used to compute the biomass CV at different scales defined in equation ( 21).

The correlation matrix R * of the system, whose elements ρ i (k) j (ℓ) are defined by:

ρ i (k) j (ℓ) = w i (k) j (ℓ) √ w i (k) i (k) √ w j (ℓ) j (ℓ) (22)

S1-8 Asymptotic resilience

In addition to the response to stochastic perturbations, we consider asymptotic resilience to measure the long term return time of the metacommunity. Asymptotic resilience is measured by the opposite of the real part of the dominant eigenvalue λ dom of Jacobian matrix J (-ℜ(λ dom )). Since the dominant eigenvalue is the eigenvalue with the largest real part and we only consider ecosystems at equilibrium (i.e. all eigenvalues have negative real parts), the lower the real part of the dominant eigenvalue, the faster the long term return time.

Moreover, we can assess the influence of each species on asymptotic resilience by comparing the absolute value of the real part of each element e i of the dominant eigenvector (E dom ). Because e i is the contribution of species i to E dom , |e i |/ n j=1 |e j | is the relative weight of species i in the dynamics of long term return to equilibrium (with n the number of populations in the metacommunity). Figure S2-2: Distribution of biomass in an isolated predator-prey system (without dispersal) depending on interaction strength a relative to self-regulation. Increasing the asymmetry of the interaction strength γ is equivalent to increasing a (m = 0.65).

S2 Complementary results

S2-1 General description of parameters

Increasing the interaction strength a relative to self-regulation decreases the biomass of prey because of the increased mortality due to predation (Figure S2-2). However, the biomass of predators follows a hump-shaped relationship: it first increases due to the increased resource consumption and then decreases because of prey overexploitation.

The effect of perturbations of populations within a community of S species can be assessed by the ratio of the mean variance of species biomass B j to the variance of perturbations σ k :

1 S S j Var(B j ) 1 S S k σ 2 k = 1 S S j Var(B j ) σ 2 i
because we only conisder one perturbation effecting species i (23)

As demonstrated by Arnoldi et al. (2019), exogenous perturbations affect more rare species, demographic perturbations evenly affect species regardless on the biomass distribution and environmental perturbations affect more abundant populations. Therefore, we consider demographic perturbations to perturb the entire community with the same intensity regardless on the biomass variations caused by varying γ (Figure S2-3). 23)) depending on the biomass of the perturbed species. Three types of perturbations with different scaling with the equilibrium biomass of the perturbed species i (B * z i ) are tested: exogenous perturbations (z = 0), demographic perturbations (z = 0.5) and environmental perturbations (z = 1). The asymmetry of interaction strength alters the biomass distribution between the two patches. Increased the interaction strength in patches #1 enhances the biomass production of predators that spill over patches from patch #1 to patch #2, therefore increasing prey biomass in patch #1 and decreasing it in patch #2 (γ = 3). Plain rectangles represent the biomass of each population, while dashed rectangles represent the same population in a metacommunity without dispersal (no spillover effect)

S2-2 Dispersal of predators and perturbation of prey

S2-2-1 Source-sink effect

Varying the asymmetry of interaction strength γ is equivalent to varying the interaction strength a in patch #1 and has the same effects on species biomass: first increasing γ increases predator biomass by increasing prey consumption, then it decreases predator biomass because of resource overexploitation (Figure S2-2). This leads to different biomass distributions in patches #1 and #2 (Figure S2-4A). Predator biomass increases with γ and is the same in both patches because their high dispersal rate balances any difference. Prey biomass is higher in patch #1 than in patch #2, and both decrease with γ, except in patch #1, where we first observe a small increase for γ < 2. This response is due to source-sink effects: the increase in prey consumption in patch #1 increases predator biomass (source) that spills over patch #2 (sink) due to dispersal (Figure S2-4B and Figure S2-5). Therefore, predator biomass is lower in patch #1 and higher in patch #2 compared to what we expect in the same food chains in isolation (i.e., without dispersal). This also prevents predators from overexploiting prey in patch #1 by spreading the increased predator biomass across the metacommunity, which explains why we do not observe a decrease in predator biomass for high values of γ, as shown in Figure S2-2. Conversely, the distribution of prey biomass across the two patches is opposite (higher in patch #1 and lower in patch #2).

S2-2-2 Conditions of coexistence

Asymmetry and dispersal lead to competition, apparent competition and source-sink dynamics that can rescue or drive local populations to extinction. Therefore, we consider limit cases in which dispersal is infinite (well mixed populations across the metacommunity) to analytically calculate biomasses at equilibrium and determine the range of values of ω and γ enabling the coexistence of all populations of each species.

We consider the total biomass of predators

B tot 2 = B (1) 2 + B (2) 2
and because the very high dispersal of predators equally distributes its biomass among the two patches, we have

B * tot 2 = 2 × B * (1)
2

. Then, we can define the system:

dB (1) 1 dt = DB (1) 1 ωg D -B (1) 1 -γma B tot 2 2 (24a) dB (2) 1 dt = DB (2) 1 g D -B (2) 1 -ma B tot 2 2 (24b) dB tot 2 dt = mD B tot 2 2 - r D -B tot 2 + εa(γB (1) 1 + B (2) 1 ) (24c)
Since r = 0, we remove it from the equations for the sake of simplicity. We define λ = εma 2 , which is the intensity of top-down control defined by Barbier and Loreau (2019). At equilibrium, we obtain:

B (1) * 1 = g D 2ω + ωλ -γλ 2 + λ(γ 2 + 1) (25a) B (2) * 1 = g D λγ 2 -ωλγ + 2 2 + λ(γ 2 + 1) (25b) B tot * 2 = 2εag(1 + ωγ) D(2 + εa 2 m(γ 2 + 1)) (25c)
Prey biomass in patch #1 B * ( 1)

1
is positive only if:

γ < ω(2 + λ) λ ----→ λ→∞ ω (26) Prey biomass in patch #2 B * (2) 1 is positive if f (γ) = λγ 2 -ωλγ + 2 > 0. f opens upwards: thus, if ω < 8/λ, f has no roots and is always positive. Otherwise, B (2) 1 is positive if: γ > λω + λ(λω 2 -8) 2λ ----→ λ→∞ ω if ω > 8 λ or (27a) γ < λω -λ(λω 2 -8) 2λ ----→ λ→∞ 0 if ω > 8 λ (27b)
Predators B 2 thrive in each patch for all values of ω and γ (Figure S2-6C). Hence, coexistence is ensured for all values of top-down control λ, asymmetry of resource supply ω and asymmetry of interaction strength γ only if γ = ω. In the main text, we always consider γ = ω (Figure S2-6D), but their independent effects are detailed in the following.

S2-2-3 Nontransitivity of correlation

To explain the correlation between prey populations, we can track the transmission of perturbations in the metacommunity. Increasing the asymmetry of interaction strength γ tends to decorrelate predator and prey dynamics within each patch (Figure S2-7A). When prey are perturbed in patch #1, the dynamics of predator and prey biomass are correlated in patch #1 and anticorrelated in patch #2 due to the bottom-up and top-down transmissions of perturbations, respectively. Although we would expect the two populations of prey to be anti-correlated according to the mechanism described by Quévreux et al. (2021) (see Figure S2-28 in the following), we actually observe a weak correlation of these two populations (Figure S2-7B). In the same way, the intermediate correlation and anti-correlation of predator and prey when prey are perturbed in patch #2 do not explain the strong anti-correlation of prey populations (Figure S2-7C). Therefore, other mechanisms are acting in our system.

S2-2-4 Complete effects of εa and ma

The stability patterns observed in Figures 2, 3 and 5 in the main text are also observed for a wide range of ecological and physiological parameters aggregated into the positive effect of prey on predators εa and the negative effect on predators on prey ma. Therefore, our results are robust and the identified mechanisms are specific to a particular combination of parameters.

The response of the asymptotic resilience to the asymmetry of interaction strength (Figure S2-11A) is not similar to the results of Rooney et al. (2006). Indeed, we do not observe minimum of resilience for γ = 1 for all combinations of εa and ma. The variations in asymptotic resilience depend on the relative contribution of each population of each species (Figure S2-11B), which is governed by the biomass distribution of each species among patches (Figure S2-8A) and the ratio of predator to prey metabolic rate ratio m (Figure S2-1B). As demonstrated by Haegeman et al. (2016) and Arnoldi et al. (2018), rare species control the long term response to perturbations of the metacommunity (i.e., the asymptotic resilience), as well as species with a slow pace of life (i.e., a slow metabolism).

S2-2-6 Effect of perturbation of predators

The perturbation of predators leads to the same response regardless of the perturbed patch because the very high dispersal of predators perfectly synchronises their population dynamics. The asymmetry of interaction strength leads to different dynamics in each patch that decreases the correlation of prey dynamics (Figure S2

S2-2-7 Effect of food chain length

Here, we consider three trophic levels to extend our results to metacommunities with longer food chain lengths. In this setup, only top predators (species 3) are able to disperse, and basal species (species 1) receive stochastic perturbations. γ also has the same value across trophic levels. We observe a similar response to the case with two trophic levels for the correlations of the dynamics of the biomass of species 2 and 3 (Figure S2-18) as well as for biomass CV . However, the response of species 1 in completely different. Therefore, the mechanisms described in the main text are only acting for the dispersing species and the species directly interacting with it. S2-18: Correlation between populations in a three trophic level food chain depending on the asymmetry of interaction strength γ , positive effect of prey on predator εa and negative effect of predator on prey ma (d 3 = 10 6 and ω = 1). Top predators disperse, and the basal species is perturbed in patch #1 or #2.

  ures S2-8 and S2-10 in the supporting information).

1 ⃝,Figure 3 :Figure 4 :

 134 Figure3: Stability at different scales depending on asymmetry of interaction strength γ when predators disperse and prey are perturbed in patch #1 or #2. A) The temporal variability in the metacommunity is assessed by the coefficient of variation (CV) of biomass at different scales: population scale (biomass CV of one species in one patch), metapopulation scale (CV of the total biomass of one species across patches) and metacommunity scale (CV of the total biomass of the entire metacommunity). B) Biomass CV at the population scale. Note that the curves for predators overlap because their high dispersal perfectly balances their biomass distribution between the two patches. C) Biomass CV at the metapopulation scale. D) Biomass CV at the metacommunity scale.

Figure 5 :

 5 Figure5: A) Time series of biomasses rescaled by their value at equilibrium after an increase of 20% in prey biomass in patch #1 (left panel) or patch #2 (right panel) for a value of interaction strength asymmetry γ = 3. B) Direct effect of prey on predator (and vice versa) depending on interaction strength asymmetry γ. Direct effects correspond to the terms of the Jacobian matrix and the dashed line emphasises the null direct effects. On the central schema, the numbers are the numeric values of the terms of the Jacobian matrix corresponding to each arrow for γ = 3.

Figure S2- 1 :

 1 FigureS2-1: Distribution of parameters and their effects on an isolated food chain. A) Biomass distribution depending on the positive effect of prey on predator ϵa and the negative effect of predator on prey ma. B) Value of the ratio of predator to prey metabolic rate m = m i+1 /m i for each combination of ϵa and ma.

Figure

  FigureS2-3: Ratio of the mean variance of species biomass to the mean variance of environmental perturbations (see equation (23)) depending on the biomass of the perturbed species. Three types of perturbations with different scaling with the equilibrium biomass of the perturbed species i (B * z i ) are tested: exogenous perturbations (z = 0), demographic perturbations (z = 0.5) and environmental perturbations (z = 1).

Figure

  Figure S2-4: A) Distribution of the biomass of each species among patches depending on the asymmetry of interaction strength γ. B) Distribution of biomasses scaled by their value in a metacommunity without dispersal (B scaled = B d2>0 /B d2=0 ).

  FigureS2-5: The asymmetry of interaction strength alters the biomass distribution between the two patches. Increased the interaction strength in patches #1 enhances the biomass production of predators that spill over patches from patch #1 to patch #2, therefore increasing prey biomass in patch #1 and decreasing it in patch #2 (γ = 3). Plain rectangles represent the biomass of each population, while dashed rectangles represent the same population in a metacommunity without dispersal (no spillover effect)

  Figure S2-6: Distribution of parameters, asymmetry of resource supply ω and asymmetry of interaction strength γ, leading to the coexistence of predator and prey in each patch. Only predators are able to disperse at an infinite rate (well-mixed predator populations). This distribution is assessed for different values of the positive effect of prey on predator ϵa and negative effect of predator on prey ma. A) Biomass of prey in patch #1 B because predator populations are well mixed). D) Coexistence of predator and prey in each patch.

Figure

  Figure S2-7: Correlation between predator and prey with each patch (εa = 1, ma = 1, d 2 = 10 6 , ω = γ). A) Correlation in each patch depending on asymmetry of interaction strength γ when prey in patch #1 (left panel) or in patch #2 (right panel) are perturbed. Labels and vertical dashed lines represent the correlation values used in panels B and C. B) Schematic representation of correlations when γ = 7. Coloured double arrows and their associated number (see also the colour scale) represent the correlations between populations. Prey are perturbed in patch #1. C) Prey are perturbed in patch #2.

  Figure S2-8: Biomass distribution of each species in each patch depending on asymmetry of interaction strength γ, positive effect of prey on predator εa and negative effect of predator on prey ma. A) Biomass distribution. B) Biomass scaled by the biomass in the metacommunity without dispersal (d 2 = 0).

Figure

  FigureS2-9: Correlation between populations depending on asymmetry in interaction strength γ when predators disperse and prey are perturbed in patch #1 or #2. Predators have a high scaled dispersal rate (d 2 = 10 6 ), which strongly couples their two populations (γ = ω).

  FigureS2-10: Biomass CV at different scales depending on asymmetry of interaction strength γ, positive effect of prey on predator εa and negative effect of predator on prey ma (d 2 = 10 6 and ω = γ). A) Biomass CV of the population of each species in each patch. B) CV of the total biomass of each species.

FigureFigure

  FigureS2-11: Linear stability depending on asymmetry of interaction strength γ, positive effect of prey on predator εa and negative effect of predator on prey ma. A) Asymptotic resilience (real part of the dominant eigenvalue of the Jacobian matrix) B) Contribution of the populations of each species to the dominant eigenvector.

Figure

  FigureS2-13: Time series of biomasses rescaled by their value at equilibrium after an increase in prey biomass by 20% in patch #1 (left panel) or patch #2 (right panel) depending on asymmetry of interaction strength γ, positive effect of prey on predator ϵa and negative effect of predator on prey ma (γ = 3, d 2 = 10 6 and ω = γ).

  -16) and stabilises predator dynamics by decreasing their biomass CV (Figure S2-17A).

Figure

  FigureS2-16: Correlation between populations depending on asymmetry of interaction strength γ, positive effect of prey on predator εa and negative effect of predator on prey ma. Predators disperse and are perturbed in patch #1 or patch #2 (d 2 = 10 6 and ω = γ).

  FigureS2-17: Biomass CV at different scales depending on asymmetry of interaction strength γ, positive effect of prey on predator εa and negative effect of predator on prey ma. Predators disperse and are perturbed in patch #1 or patch #2 (d 2 = 10 6 and ω = γ). A) CV of the total biomass of each species. B) CV of the total biomass of the metacommunity.

Figure

  FigureS2-18: Correlation between populations in a three trophic level food chain depending on the asymmetry of interaction strength γ , positive effect of prey on predator εa and negative effect of predator on prey ma (d 3 = 10 6 and ω = 1). Top predators disperse, and the basal species is perturbed in patch #1 or #2.

FigureFigure

  FigureS2-20: Distribution of parameters, asymmetry of resource supply ω and asymmetry of interaction strength γ, leading to the coexistence of predator and prey in each patch. Only prey are able to disperse at an infinite rate (well mixed prey populations). This distribution is assessed for different values of the positive effect of prey on predator εa and the negative effect of predator on prey ma. The product of εa and ma is the strength of top-down control λ (λ = εa 2 m, seeBarbier and Loreau (2019)). A) Biomass of prey in patches #1 and #2 (B

  FigureS2-23: Biomass CV at different scales depending on asymmetry of interaction strength γ, positive effect of prey on predator εa and negative effect of predator on prey ma (d 1 = 10 6 and ω = γ). A) Biomass CV of the population of each species in each patch. B) CV of the total biomass of each species.

  Figure S2-24: Time series after pulse perturbation of predators in patch #1 or in patch #2. Biomasses are scaled by their value at equilibrium, and dispersal is high (γ = 3 and d = 10 6 ).

  for an extended review and more references.

	Predator	Prey	γ	Reference
	American marten	Vole species	1.6	Andruskiw et al., 2008
	Wolf	Moose	14-100	Gervasi et al., 2013
	Wolf	Roe deer	2.5-8	Gervasi et al., 2013
	Wolf	Elk	10	Kauffman et al., 2007
	Savannah predators Savannah ungulates 1.5-4.5	Thaker et al., 2011
	Artificial gecko	Australian predators	2.8	Hansen et al., 2019
	Lynx	Roe deer	2	Gehr et al., 2020
	Puma	Vicuña	1.6	Donadio and Buskirk, 2016

  .

	Parameter	Interpretation	Value
	σ i	standard deviation of stochastic noise 10 -3
	g	net growth rate of prey	1
	r	death rate of predators	0
	D	self-regulation	1
	ϵ	conversion efficiency	0.65
	m	predator/prey metabolic rate ratio	0.65
	a	attack rate	1.54,
	εa	positive effect of prey on predators	1
	ma	negative effect of predators on prey	1
	d 2	scaled dispersal rate of predators	10 6
	ω	asymmetry of resource supply	[1,10]
	γ	asymmetry of interaction strength	[1,10]

2 1 perturbation of prey in patch #1 prey in patch #2

Acknowledgements

This work was supported by the TULIP Laboratory of Excellence (ANR-10-LABX-41) and by the BIOSTASES Advanced Grant, funded by the European Research Council under the European Union's Horizon 2020 research and innovation programme (666971).

This paper has been peer-reviewed and recommended by PCI Ecology and Werner Ulrich (10.24072/pci.ecology.100512). We thank the three reviewers Phillip P.A. Staniczenko, Ludek Berec and Diogo Provete for their constructive comments.

Data accessibility

The R codes to reproduce the results and the figures are available on GitHub (https://github.com/PierreQuevreux/model_ metacommunity_spatial_heterogeneity).

A B

Figure S2-19: Biomass CV in a three trophic level food chain at different scales depending on asymmetry of interaction strength γ, positive effect of prey on predator εa and negative effect of predator on prey ma (d 3 = 10 6 and ω = 1). A) Biomass CV of the population of each species in each patch. B) CV of the total biomass of each species.

S2-3 Dispersal of prey and perturbation of predators

In this section, we consider a setup mirroring the metacommunity model described in the main text. Here, only prey are able to disperse at a very high rate (d 1 = 10 6 ), and predators receive stochastic perturbations. We also set γ = ω to be consistent with the results in the main text. In the following, we find the same responses to the asymmetry of interaction strength γ, which demonstrates that the mechanisms described in the main text are not conditioned by the trophic position of the dispersing species.

S2-3-1 Conditions of coexistence

When only prey are able to disperse, all populations of each species have positive biomasses for all values of ω, γ, εa and ma (Figure S2-20). Here, we consider the same metacommunity as in the main text (see Figure 1 in the main text), but prey receive spatially correlated environmental perturbations. Environmental perturbations correspond to the synchronous response of all individuals of the same population to an environmental factor (e.g., drought), and they scale with equilibrium biomass B * i (see the supporting information of Quévreux et al. (2021) for the demonstration). In our metacommunity, we also consider that environmental perturbations are spatially correlated since it is reasonable to assume that different populations of the same species will respond in a similar way to environmental perturbations.

S2-4 Correlated environmental perturbations

The effect of a perturbation on a population within a community can be assessed by the ratio of the mean variance of species biomass j to the variance of the perturbation i by equation ( 23). As demonstrated by Arnoldi et al. (2019), environmental perturbations affect abundant populations the most, which is the prey population in the fast patch in our case . Therefore, we can approximate the effect of environmental perturbations by the effect of the perturbation of prey in the fast patch. Perturbing a single population with demographic or environmental perturbations leads to exactly the same qualitative results (Figure S2 Increasing the correlation of perturbations increases the correlation of the dynamics of prey populations (Figure S2-27A) because of the Moran effect [START_REF] Moran | The statistical analysis of the Canadian Lynx cycle[END_REF]. The increase in synchrony explains the increase in the biomass CV observed at each scale for all species (Figure S2-27B-D), except for prey in the slow patch (Figure S2-27B). The Moran effect is particularly strong at low asymmetry (γ < 4), but once asymmetry is high enough, two mechanisms disrupt the Moran effect. First, when asymmetry is high, the dynamics in each patch become so different that correlated perturbations are not able to generate similar responses. Second, because of the discrepancy in the distribution of prey biomass among the two patches, environmental perturbations mostly affected prey in the fast patch . Therefore, with increasing asymmetry of interaction strength γ, the response of the metacommunity to correlated environmental perturbations converges towards the response of a metacommunity in which only prey in the fast patch are perturbed. 2021), who considered a two patch predator-prey metacommunity with passive dispersal. In the setup presented in A), prey are perturbed in patch #1 and only predators are able to disperse. Thus, perturbations have a bottom-up transmission in patch #1 (i.e. transmission from lower to upper trophic levels). This leads to the temporal correlation of the biomass dynamics of predators and prey in patch #1 showed in B)(1) because if a perturbation increases the biomass of prey, it also increases the biomass of predators due to the vertical transfer of biomass. The passive dispersal of predators transmits the perturbations and spatially correlate their populations as shown in B)(2). Then, perturbations have a top-down transmission in patch #2 (i.e. transmission from upper to lower trophic levels). This leads to the temporal anticorrelation (negative coefficient of correlation) of the biomass dynamics of predators and prey in patch #2 showed in B)(3) because if a perturbation increases the biomass of predators, it decreases the biomass of prey due to the negative effect of predators on prey. Eventually, prey populations are spatially anticorrelated, as shown in B)(4). Hence, by knowing which species is perturbed, which species disperses and how perturbations propagate within a food chain, Quévreux et al. (2021) were able to explain the spatial synchrony of the various populations of a metacommunity, summariesed by the correlation matrix in C).