Introduction

Anthropogenic habitat degradation and loss are strong negative drivers of biodiversity on local and global scales (Butchart et al. 2010; Pereira et al. 2010; Pimm et al. 2014). While habitat loss has a clear cause–effect relationship with declining diversity induced by, e.g., lack of resources, habitat size restrictions or increased mortality (Brooks et al. 2002; Duraiappah et al. 2005), the effect of other modifications of the landscape such as fragmentation is still intensely debated (Hanski 2015; Fahrig 2017; Fletcher et al. 2018; Fahrig et al. 2019). Following Fahrig (2003), habitat fragmentation comprises three main components: the number of patches, patch isolation and patch size, but excludes habitat loss. Their respective effects are more difficult to assess because they are usually weaker than the effects of habitat loss (Fahrig 2003) and often confounded with the latter (Didham et al. 2012).

In metacommunities, patch isolation determines to which extent individuals can disperse through the landscape and thereby contribute to the regional distribution and persistence of species. Empirical and experimental studies report, however, conflicting results of patch isolation at different spatial scales: Negative effects on regional diversity have been attributed to the prevention of rescue effects (Levins 1969; Gotelli 1991), but also positive effects on local diversity have been recorded (Fahrig 2017). On the local scale dispersal can also alter biotic interactions among species directly, emphasizing the interplay between local and regional dynamics in metacommunities (Walting and Donnelly 2006). Recent modeling approaches on metacommunities try to integrate more details of local and regional aspects regarding landscape attributes and species interactions, but mainly focus on species persistence and diversity (Pillai et al. 2011; Ryser et al. 2019) and ignore effects of dispersal on local population dynamics and its relevance for stability (LeCraw et al. 2014).

A major concern of models that include explicit population dynamics are mechanisms that synchronize population cycles between habitat patches. Such synchronous oscillations destabilize metapopulations by amplifying the amplitude of oscillations in their regional abundances and increasing the extinction risk of species in entire regions due to correlated local extinction events. Conversely, asynchronous oscillations can promote regional persistence and stability through rescue effects (Levins 1969; Blasius et al. 1999) or the portfolio effect (Schindler et al. 2015; Thorson et al. 2018). These models, which are often limited to either a small number of patches or to regular, rectangular lattices (Briggs and Hoopes 2004), have established that the synchronicity of population oscillations between patches generally increases with dispersal rate (Sherratt et al. 2000; Jansen 2001). Other factors affecting synchronicity are adaptive dispersal (Abrams 2007; Abrams and Ruokolainen 2011), inter- and intraspecific density dependence of dispersal rates (Hauzy et al. 2010), and costliness or distance dependence of dispersal (Koelle and Vandermeer 2005). In larger networks of habitat patches, an irregular network structure favors asynchronous dynamics (Holland and Hastings 2008), but high dispersal rates again lead to synchronous oscillations that are detrimental for species persistence (Plitzko and Drossel 2015). At larger effective distance between patches, dispersal between them is limited (Koelle and Vandermeer 2005; Fletcher et al. 2016), linking the results regarding synchronization of population oscillations to research on the effect of patch isolation. Indeed, it has been shown that synchronization among natural populations declines with increasing distance between them (Ranta et al. 1995).

While synchronization is often linked to dispersal rate and thereby implicitly to landscape properties like patch isolation, it can also be directly affected by correlated environmental fluctuations (Moran 1953; Ranta et al. 1995; Koenig 1999; Kahilainen et al. 2018). These fluctuations can affect demographic rates of the species via changing environmental conditions (like ambient temperature or resource availability), but they can also directly influence the availability of patches as habitable areas. As an example for the latter, a landscape in which both a temporally variable environment and a pronounced spatial structure strongly affect ecological communities is kettle holes in formerly glaciated regions (Kalettka and Rudat 2006). These small ponds are typically formed in large clusters, and seasonal changes of temperature and precipitation cause some of them to be only temporally filled with water. The local aquatic communities of these temporary ponds thus periodically become completely extinct, and recolonization through dispersing species from permanent ponds is a key element to reestablish the communities (De Meester et al. 2005). As the recolonization happens in a temporally correlated manner at the beginning of the wet season, a synchronizing effect on the population dynamics can be expected. However, this is again contingent on the spatial structure of the landscape, as lower dispersal rates due to higher mean patch isolation can impede the recolonization process.

So far, the interaction between these drivers of synchronization and population variability in general remains largely unexplored (but see (Gouhier et al. 2010)), despite the fact that anthropogenic activity continues to increase both habitat degradation and environmental variability. In order to fill this gap, we examine the dynamics of a meta-food chain in large, spatially explicit networks of habitat patches and analyze its stability with respect to the mean patch isolation of the landscape and environmental disturbances that periodically render a subset of the patches uninhabitable. We chose a food chain as model system because it has, on the one hand, a simple and tractable structure that, on the other hand, already allows for indirect effects mediated by feeding interactions on different trophic levels. In order to obtain a complete picture of the effects of patch isolation and periodic environmental disturbances on the extent and synchronicity of population oscillations in food chains, we analyze two parameterizations of the food chain that correspond to contrasting oscillation patterns. These patterns are characterized either by a relatively even distribution of biomass along the food chain (weak trophic cascade) or by marked differences among the species (strong trophic cascade), both of which are common in natural ecosystems (Estes and Duggins 1995; Carter and Rypstra 1995).

Our model setup explicitly addresses one aspect of fragmentation, namely patch isolation, while keeping other potentially confounding drivers such as the total amount of habitat or the number of patches constant. We consider both static landscapes, where all patches are constantly available as habitats, and dynamic landscapes, where periodic environmental disturbances regularly render some of the patches uninhabitable. The stability of the dynamics of the metacommunity is evaluated within the framework of Wang and Loreau (2014) that divides population variability into an \(\alpha\)-, \(\beta\)-, and \(\gamma\)-component (similar to the classical diversity indices by Whittaker (1972)): \(\alpha\)-variability is the average coefficient of variation of a species’ local abundances, \(\gamma\)-variability is the coefficient of variation of the regional (metapopulation) abundance, and \(\beta\)-variability quantifies differences in oscillations between patches, i.e., how synchronously the local populations oscillate. Generally, it is assumed that higher dispersal rates synchronize population dynamics (e.g., Gouhier et al. (2010)). When mean patch isolation increases, mortality during dispersal increases, too. We expect that this decreases net dispersal flows and thus also decreases synchrony of population dynamics among patches (i.e., increases \(\beta\)-variability). This may, however, be counteracted by (synchronous) periodic disturbances of patch availability. Furthermore, we expect local (\(\alpha\)-) variability to decrease, as increasing mortality allows less biomass to flow up the food chain, thus weakening (and thereby stabilizing) the trophic interactions (Rip and McCann 2011). If the local population oscillations indeed become less synchronous, this will also decrease regional (\(\gamma\)-) variability as habitats become more isolated.

Methods

The model comprises a tri-trophic food chain including an autotroph (A), a consumer (C) and a predator (P) species. As basis for the growth of the autotroph, a dynamic resource (R) serves as essential energy source and can be seen as a universal nutrient. This food chain is extended to a metacommunity by placing copies of it on habitat patches that are randomly distributed in space and connected via species-specific dispersal links (Fig.  1). Where applicable, the individual parameters are derived from empirical data, largely from invertebrate communities.

Fig. 1
figure 1

a) simplified example of a spatial network of habitat patches. Dashed lines of different grey tones indicate dispersal links of the respective species. The resource does not disperse between patches. b) local food chain on each patch comprising three trophic levels (autotrophs, A, consumers, C, and predators, P) plus a dynamic resource, R

Full size image

Trophic interactions

We first describe only the trophic interactions between the populations on a single patch and disregard dispersal. The local dynamics of the food chain follow a generalization of the bioenergetics approach (Yodzis and Innes 1992; Brose et al. 2006), supplemented with an equation for the resource. Adapted from chemostat dynamics, the rate of change of the resource density R is expressed as

$$\begin{aligned} \frac{dR}{dt} = D\cdot \left( R_{0}-R\right) -G_{AR} A \end{aligned}$$
(1)

with the resource turnover rate D and the supply concentration \(R_0\). Uptake of resources by the autotroph A is described by a Monod function \(G_{AR}=r\frac{R}{K+R}\) with maximum uptake rate r and half saturation constant K. The rates of change in biomass density for each species (A, C and P) are expressed by

$$\begin{aligned} \frac{dA}{dt}= & {} G_{AR} A-F_{CA}C-x_{A}A=g_A A\nonumber \\ \frac{dC}{dt}= & {} e_{C}F_{CA}C-F_{PC}P-x_C C=g_C C \nonumber \\ \frac{dP}{dt}= & {} e_{P}F_{PC}P-x_{P}P=g_P P \end{aligned}$$
(2)

where the first terms in all three equations represent growth due to consumption, the last terms denote metabolic losses, and the middle terms in the equations for the autotroph and the consumer describe mortality through predation. The terms are summarized by the net per capita growth rates \(g_i\) (\(i=A,C,P\)). The parameters \(e_i\) and \(x_i\) are assimilation efficiencies and per capita respiration rates, respectively. The per capita feeding rate of species i on species j is described by a Beddington–DeAngelis functional response (DeAngelis et al. 1975; Beddington 1975):

$$\begin{aligned} F_{ij} = \frac{1}{m_{i}} \frac{a_{ij}B_{j}}{1+a_{ij}h_{ij}B_{j}+c_{i}B_{i}} \end{aligned}$$
(3)

with the attack rate \(a_{ij}\), the handling time \(h_{ij}\), the interference coefficient \(c_i\), and \(B_{i}\) and \(B_j\) as placeholders for the respective consumer’s or resource’s biomass density. Since the model is formulated in terms of biomass densities (as opposed to population densities), the functional response is scaled with \(\frac{1}{m_i}\), the inverse of the respective consumer’s body mass (Heckmann et al. 2012).

The parameters of the trophic dynamics scale allometrically with the body mass of the species. Mass-specific maximum growth rate and respiration rates are assumed to decrease with a negative quarter-power law with body mass, i.e., \(r=r_0m_i^{-0.25}\) and \(x_i=x_{0,i}m_i^{-0.25}\) (Yodzis and Innes 1992; Brose et al. 2006). Following Rall et al. (2012), handling times depend on the body masses of both consumer and resource with \(h_{ij}=h_0m_i^{-0.48}m_j^{-0.66}\). The same is true for the attack rates, but since these parameters were used to differentiate the contrasting states of top-down control, fixed values were used here (c.f. Table 1) that nevertheless obey the general trends found in Rall et al. (2012). Body masses increase by a factor of 100 per trophic level, a value commonly found in invertebrate communities and known to have a stabilizing effect on population dynamics (Brose et al. 2006; Brose et al. 2006b). Freedom of choosing an appropriate set of units allows us to set the body mass of the autotroph to \(m_A=1\). In general, the model is parameterized such that the population dynamics of all species are oscillatory when dispersal is not accounted for (Table 1, Fig. 2).

Habitat network and dispersal

We use the same rules for modeling spatial interactions as in Ryser et al. (2019). Dispersal is considered for the autotroph, consumer and predator species in the model. The spatial setting is implemented as a random geometric graph (RGG) (Penrose 2003), where each node of the spatial network represents a habitat patch for a local community (Urban and Keitt 2001). The (x, y)-coordinates of each patch were drawn at random from a bivariate uniform distribution over the interval \(\left[ 0:1\right] \times \left[ 0:1\right]\). Dispersal links between the patches connect the local populations, enabling exchange of biomass between patches and thereby forming a meta-food chain (Fig.  1).

Each species perceives its individual dispersal network depending on its body mass \(m_{i}\). A dispersal link for species i exists between two patches k and l only if the distance between them is less than the species-specific maximum dispersal distance

$$\begin{aligned} D_{max,i}=D_{0}m_{i}^{\epsilon }\,. \end{aligned}$$
(4)

The exponent \(\epsilon\) is set to a positive value to account for increased mobility and thus improved dispersal abilities of species with a larger body mass (Hein et al. 2012; Peters 1983).

Dispersal itself is at least for animal species often an active process resulting in metabolic costs and potentially involving a higher risk of predation. To account for these costs (dispersal mortality), we assume that dispersal success \(S_{i,lk}\) (i.e., the fraction of individuals not dying during dispersal) of species i, when moving between patches l and k, decreases linearly with the distance between the patches:

$$\begin{aligned} S_{i,lk}=\max \left( 1-d_{i,lk}, 0\right) \,, \end{aligned}$$
(5)

where \(d_{i,lk}=\frac{d_{lk}}{D_{max,i}}\) is the distance between the patches relative to the maximum dispersal distance of species i. For passively dispersing plants, distance-depending costs can be caused by a decreasing probability of propagules finding by chance a suitable patch that is further away.

The fraction of individuals emigrating from a source patch k that move toward a target patch l is calculated using the weight function

$$\begin{aligned} W_{i,lk} = \frac{1-d_{i,lk}}{\sum _{p}(1-d_{i,pk})}\,, \end{aligned}$$
(6)

where the sum in the denominator is taken over all potential target patches p that are within the maximum dispersal range of species i on patch k (i.e., those with \(d_{pk}<D_{max,i}\)). This weight function makes dispersal links between nearby patches stronger, implying that a larger proportion of emigrating biomass arrives there, than those between patches that are further apart. Note that while specific distances \(d_{i,lk}\) and success terms \(S_{i,lk}\) are symmetric for all pairs of patches, the weight function is not (i.e., \(W_{i,lk}\ne W_{i,kl}\)).

In general, the process of dispersal can be described as an exchange of biomass between habitat patches that is affecting the population dynamics of species i on patch l via emigration (\(E_{i,l}\)) from this patch and immigration (\(I_{i,l}\)) into the patch. The full population dynamics of species i on patch l, comprising both local, trophic dynamics, Eqs. (2), and dispersal dynamics, can thus be written as

$$\begin{aligned} \frac{dB_{i,l}}{dt}=g_{i,l}B_{i,l} - E_{i,l} + I_{i,l}\,. \end{aligned}$$
(7)

Emigration is a complex process in nature possibly involving different environmental cues and species properties. Here, we assume an adaptive emigration rate that depends on the net per capita growth rate \(g_{i,l}\) of species i on patch l, reflecting its current situation in this habitat. If a species’ net growth is positive, there is little need for dispersal and emigration will be low. However, if the local environmental conditions deteriorate, e.g., due to low resource availability or high predation pressure, the emigration rate increases. This is captured by the following function:

$$\begin{aligned} E_{i,l} = \frac{\mu _{i}B_{i,l}}{1+e^{b\left( g_{i,l}+x_{i}\right) }}\,. \end{aligned}$$
(8)

The parameter \(\mu _{i}=\mu _0 x_i\) determines the maximum per capita emigration rate and b determines how sensitively the emigration rate depends on the net growth rate (i.e., how quickly it drops when \(g_{i,l}\) increases). Finally, immigration of species i into patch l depends on the amount of emigration from all neighboring patches k as well as on the specific dispersal network, encoded in the success and weight functions \(S_{i,lk}\) and \(W_{i,lk}\), according to

$$\begin{aligned} I_{i,l} = \sum _{k} S_{i,lk} W_{i,lk} E_{i,k}\,. \end{aligned}$$
(9)

The parameters defining the dispersal dynamics are also summarized in Table 1.

Table 1 Standard parameter set used in the model
Full size table

Simulation setup

Static and dynamic landscapes

The baseline simulations are carried out using static landscapes, i.e., with RGG networks of \(Z=30\) habitat patches as described above, where all patches and dispersal links are permanently available. However, since the environmental conditions in nature are rarely completely constant, we also study dynamic landscapes in which a fraction \(\sigma\) of the patches becomes periodically unavailable as a habitat. This process is called “blinking” and has a period length \(\lambda =6000\). This period length encompasses several hundred generation times of the autotroph, thereby providing sufficient time for the food chain to recover between blinking events. Blinking patches are turned on and off synchronously and change their state every \(\frac{\lambda }{2}\) time units. When the blinking patches are turned off, the local food chains go extinct immediately. Furthermore, the dispersal network can be disrupted because these patches cannot be used as stepping stones for dispersal between patches that are too far apart for a direct dispersal link.

Patch isolation

To capture the effects of varying mean patch isolation, the intercept of the maximum dispersal distance, \(D_0\), (Eq. (4)) is varied systematically between 0.06 and 0.5. This creates habitat networks that range from mostly isolated patches to systems where the predator can move in a single step between any two patches. The spatial network is quantified by the mean patch isolation of the predator’s dispersal network,

$$\begin{aligned} I_{RGG,P} = 1 - \frac{L_{P}}{\frac{1}{2}Z\cdot \left( Z-1\right) }\,, \end{aligned}$$
(10)

with \(L_P\) the number of undirected dispersal links of the predator and Z the number of habitat patches. Note that using the isolation of the dispersal network of any of the other species to define the mean patch isolation of the landscape would only rescale the x-axis of the results (Fig. 3), but not change them qualitatively.

Ecosystem stability

We evaluated ecosystem stability according to Wang and Loreau (2014) as \(\alpha\)-, \(\beta\)-, and \(\gamma\)-variability of autotroph, consumer and predator. For the mean local or \(\alpha\)-variability of a species, the coefficients of variation (CV, \(\frac{\text {standard deviation}}{\text {mean}}\)) of its local biomass densities on all patches are calculated and then averaged across patches (weighted with the respective local mean biomass density), while for the \(\gamma\)-variability (variability of the metapopulation) the CV of the total biomass density (sum over all patches) is evaluated. Similar to the \(\alpha\)-, \(\beta\)-, and \(\gamma\)-diversity indices (Whittaker 1972), \(\beta\)-variability measures differences between the patches and can thus be used to determine how synchronously local biomass densities on the different patches oscillate. It is here defined as \(\beta =\frac{\alpha }{\gamma }\). In contrast to the diversity indices, however, variability decreases with an increase in spatial scale, i.e., \(\gamma \le \alpha\) and thus \(\beta \ge 1\). Spatially synchronous oscillations result in a low \(\beta\)-variability and a \(\gamma\)-variability that approaches the value of the \(\alpha\)-variability. Perfect synchronicity is obtained at \(\beta =1\). The variability measures of a species do not change if it is permanently extinct on one or several patches. An intuitive example of two species, one with synchronous and one with asynchronous oscillations, is provided in the Online Resource (Fig. S1).

Numerical simulations

We simulated food chains that were parameterized to exhibit either a strong or a weak trophic cascade, corresponding to a very uneven or a relatively even distribution of biomass along the food chain, respectively. The weak trophic cascade was generated by relatively low attack rates of the consumer and predator species (\(a_{CA}=105\), \(a_{PC}=450\), Fig. 2a), while for the strong trophic cascade much higher attack rates were chosen (\(a_{CA}=170\), \(a_{PC}=10 000\), Fig. 2b). The spatial networks were either static (all patches permanently available as habitats) or dynamic (\(30\%\) of the patches periodically becoming unavailable as habitats). The mean patch isolation was constant for each individual simulation run, but was gradually varied between simulations by decreasing \(D_0\) from 0.5 to 0.06 in steps of 0.01. Simulations were carried out with a full-factorial design and 30 replicates for each combination of parameters, resulting in a total of 5400 simulation runs. Replicates differed in the randomly chosen positions of 30 patches that formed the spatial networks. Time series were simulated for 90 000 time units and split in three sections of equal length. During the first section, the systems settled on the attractor and from the second section, mean biomass densities were calculated. These mean biomass densities were then used to calculate the variability coefficients from the third section of the time series. During the simulations, a species was considered extinct on a given patch if its local biomass density fell below \(10^{-20}\). Global extinction of a species from the entire meta-food chain was never observed. Numerical simulations of the ODE model were performed in C (source code adopted from (Schneider et al. 2016)) using the SUNDIALS CVODE solver (Hindmarsh et al. 2005) with absolute and relative error tolerances of \(10^{-10}\). Output data were analyzed using Python 2.7.11, 3.6 and several Python packages, in particular NumPy and Matplotlib (Oliphant 2015; Van der Walt et al. 2011; Hunter 2007).

Results

Food chain dynamics without dispersal

To capture how different parameterizations of trophic interactions affect the metacommunity dynamics, we analyzed two contrasting trophic cascades in the food chain that were created by assuming either low or high attack rates. The first type, called weak trophic cascade, is characterized by a weak predation pressure of the predator, a relatively even distribution of biomass along the food chain and a high oscillation frequency (note the different scales of the x-axes of the two panels in Fig. 2). The strong trophic cascade is, in contrast, characterized by a very uneven distribution of biomass with a strong dominance of the autotroph (caused by the suppression of the consumer by the predator), a much lower oscillation frequency and much more drastic population cycles that drive both the predator and the consumer biomass densities repeatedly to very low values. The difference between the predator attack rates in the two cases had to be this pronounced as for intermediate values, the food chain is stable and the analysis of (meta-)population variabilities is not possible (Online Resource, Fig. S2).

Fig. 2
figure 2

Time series of the dynamics for the weak A and strong B trophic cascades on a single patch without dispersal dynamics. In case A, \(a_{CA}=105\) and \(a_{PC}=450\); in case B, \(a_{CA}=170\) and \(a_{PC}=10 000\). All other parameters are listed in Table 1. Note the different scales of x- and y-axes in the two panels

Full size image

Metacommunity dynamics

We evaluated the two different landscape scenarios (static vs. dynamic) for both the weak and strong trophic cascade over a gradient of the mean patch isolation. All scenarios are evaluated with respect to local (\(\alpha\)-variability), between patch (\(\beta\)-variability) and metapopulation dynamics (\(\gamma\)-variability). The observed trends in population variabilities on the different spatial scales were always the same for all trophic levels. We therefore only show results for the predator species. Results for the autotroph and consumer species are in the Online Resource (Figs. S3 and S4).

Local dynamics: \(\alpha\)-variability

In contrast to our expectations, increasing mean patch isolation amplifies biomass oscillations in static landscapes (increasing \(\alpha\)-variability, Fig. 3a,b). This trend is particularly pronounced in the strong trophic cascade from intermediate mean patch isolation (where many systems even settle on a stable fixed point) to high mean patch isolation (Fig. 3b). Because \(\alpha\)-variability has nonzero values at low mean patch isolation, the overall pattern is u-shaped. In the weak trophic cascade, \(\alpha\)-variability monotonously increases with mean patch isolation. In dynamic landscapes, \(\alpha\)-variability is higher than in static landscapes, but its main trends with mean patch isolation are significantly weaker than in static landscapes (cf. also Table 2).

Fig. 3
figure 3

Local (\(\alpha\)-variability, top row), between patch (\(\beta\)-variability, middle row) and metapopulation dynamics (\(\gamma\)-variability, bottom row) of the predator for the weak (left column) and the strong trophic cascade (right column). Light gray data points and dashed trend lines (second order fit) indicate static landscapes, and dark gray data points and solid trend lines indicate dynamic landscapes. Each data point represents the result of one simulation run with a unique spatial network of habitat patches. All data points where the variability is below \(10^{-6}\) are set to \(10^{-6}\) as differences between them provide no meaningful information that close to the fixed point

Full size image

Synchronization of patches: \(\beta\)-variability

On the regional scale, we evaluated to what extent the biomass dynamics between habitat patches synchronized (Fig. 3c,d). In line with our expectations, there is in most cases a clear trend toward decreased synchronization (increased \(\beta\)-variability, c.f. also Table 2) of the dynamics as mean patch isolation increases. The apparent limitation of synchronization in dynamic landscapes (minimal \(\beta\)-variability \(\approx 2\) for both weak and strong trophic cascades) is only a numerical effect due to the difference between constant and blinking patches.

Only the weak trophic cascade in static landscapes deviates from the general trend: Not only the \(\beta\)-variability is higher than in the other cases, but it also appears to decrease from low to intermediate mean patch isolation and only slightly increases at high mean patch isolation. The initial decrease is due to a separate cloud of data points with very high \(\beta\)-variabilities, which emerges for \(I_{RGG,P}\lessapprox 0.4\). This suggests that in this part of the parameter space a second attractor with even less synchronization between the patches exists. The bistability of the system is indeed confirmed by dedicated simulations using spatial networks with fixed coordinates of the patches (cf. Online Resource, Fig. S5)

Table 2 Summary of the trends of \(\alpha\)-, \(\beta\)- and \(\gamma\)-variability with increasing mean patch isolation for the weak (WTC) or strong (STC) trophic cascade in static or dynamic landscapes
Full size table

Metapopulation : \(\gamma\)-variability

For both the weak and strong trophic cascades, we find a relatively constant total biomass of the metapopulation (\(\gamma\)-variability \(<10^{-1}\), Fig. 3e,f). As expected, \(\gamma\)-variability is higher in dynamic landscapes than in static ones. Since local biomass oscillations are often highly synchronized, the trends in the metapopulation dynamics largely follow those already observed in the local dynamics (cf. also Table 2). As with the \(\beta\)-variability of the weak trophic cascade in static landscapes, at low mean patch isolation (\(I_{RGG}\lessapprox 0.4\)) a small cloud of data points appears to be separated from the rest, which have a low \(\gamma\)-variability. Again, these data points can be attributed to an alternative attractor with less synchronized dynamics and correspondingly a lower \(\gamma\)-variability.

Discussion

The impact of habitat fragmentation on biodiversity and community dynamics is a subject of ongoing debate (Fahrig et al. 2019; Fletcher et al. 2018). Here, we evaluated the effect of mean patch isolation as one aspect of fragmentation on the population dynamics of two contrasting states of a meta-food chain in static and dynamic landscapes. Most intriguingly, we found that both local (\(\alpha\)-) and metacommunity (\(\gamma\)-) variability increased with increasing mean patch isolation, despite the fact that synchronization among patches mostly decreased (\(\beta\)-variability increased) along the same gradient. Periodic environmental disturbances that rendered some patches regularly uninhabitable in dynamic landscapes weakened these trends, but at the prize of overall higher levels of \(\alpha\)- and \(\gamma\)-variability.

Interactions between dispersal and local interactions drive the dynamics in static landscapes

Higher effective dispersal rates at low patch isolation have been shown to synchronize the dynamics of metacommunities (Gouhier et al. 2010), but our results suggest that the extent of this effect may depend on the local interactions between the populations. While our results largely confirm the negative correlation between mean patch isolation (and thus, by proxy, effective dispersal rate) and synchronization, we also observe a significant deviation from this trend in the weak trophic cascade at low mean patch isolation. There, an alternative attractor with very asynchronous population oscillations (high \(\beta\)-variability) emerges. However, \(\alpha\)-variability is also relatively low on this attractor, which may explain the lack of synchronization: When the local populations do not oscillate much, their emigration rates are also almost constant over time, and there is consequently little potential for affecting the population oscillations on neighboring patches. This highlights the importance of details of the local interactions between species (in this case low attack rates in the weak trophic cascade that limit \(\alpha\)-variability) for collective phenomena like synchronization. Other theoretical studies also indicate a relevance of local interactions for the synchronization of population dynamics. Koelle and Vandermeer (2005) show, for example, opposing trends of synchronization between species in a food chain, which are due to an interaction between dispersal patterns and trophic interactions. Moreover, empirical studies provide evidence that dispersal may even alter biotic interactions between species directly (Walting and Donnelly 2006), further underlining the importance of local species interactions for our understanding of metapopulation dynamics.

Indirect effects of local trophic interactions also explain why our initial hypothesis, regarding decreasing \(\alpha\)-variability at increasing mean patch isolation, turned out to be incorrect in the weak trophic cascade. The hypothesis was based on the “principle of energy flux” (Rip and McCann 2011), according to which an increasing (dispersal) mortality at higher mean patch isolation should weaken and consequently stabilize the trophic interactions along the food chain (and thus decrease \(\alpha\)-variability). In contrast to this prediction, high dispersal mortality does not generally result in a lower \(\alpha\)- or \(\gamma\)-variability in our model. We attribute this counter-intuitive trend to an indirect effect of dispersal mortality: Despite their superior dispersal abilities, higher trophic levels often suffer most from mean patch isolation because they are energetically more limited than the species on lower trophic levels (Ryser et al. 2019). In fact, we also find that the higher the mean patch isolation, the lower the mean biomass of the predator (see Online Resource, Fig. S6). This decreases the per-capita predation mortality of the consumer, which more than compensates for the increase in the consumer’s dispersal mortality. In line with the principle of energy flux, this destabilizes the consumer–autotroph interaction. At high mean patch isolation, the \(\alpha\)-variability of the predator thus increases because the dynamics of the predator is driven by the increasingly unstable consumer–autotroph interaction.

This apparent mismatch between increasing \(\beta\)-variability (more asynchronous dynamics) and simultaneously increasing \(\gamma\)-variability at high mean patch isolation has also implications for the so-called portfolio effect (Schindler et al. 2015), which is often considered in more applied contexts. Specifically, the spatial portfolio effect (Thorson et al. 2018) measures how much \(\gamma\)-variability is reduced relative to its theoretical maximum (here given by \(\gamma\)-variability = \(\alpha\)-variability) due to asynchronous oscillations among different spatial locations. While we do observe such a reduction of \(\gamma\)-variability relative to \(\alpha\)-variability when mean patch isolation increases, the indirect effect of dispersal mortality discussed above still leads to an increase in \(\gamma\)-variability in absolute terms. This underlines that assessing factors that affect the synchronization of population dynamics across space is not always sufficient to understand the variability of a population on the regional scale.

Bistability in the weak tropic cascade

In static landscapes, the weak trophic cascade is bistable for low-to-medium mean patch isolation. In this parameter range, in addition to the attractor with intermediate synchronicity, which exists for the entire range of mean patch isolation, a second attractor with very asynchronous dynamics between the patches exists. Interestingly, the bistability concerns only the synchronicity of the dynamics (and consequently the \(\gamma\)-variability). Local (\(\alpha\)-) variability is not affected by whether the populations on different patches cycle more or less in synchrony (Fig. 3a).

Such bistability is relevant because it implies hysteresis (Scheffer et al. 1993): A small change in environmental conditions can drive the system away from one attractor, but for the system to return to it, a much larger change of the environmental conditions in the opposite direction will be necessary. This is particularly concerning here: The second attractor, which may be regarded as more desirable due to its lower metapopulation variability, loses its stability when the mean patch isolation increases beyond a certain threshold. However, the system may never return to it even when environmental conditions improve again, because the primary attractor never loses its stability.

A possible explanation for the occurrence of the alternative synchronization patterns we observe is the way the dispersal rate is modeled. Specifically, that the rate at which individuals emigrate from a given patch depends on the net growth rate they experience there. Emigration can thus be driven by a lack of resources (in which case emigration helps ending the unfavorable growth conditions and is thus self-limiting) or by an exceedingly high predation rate (in which case emigration actually intensifies the per-capita predation rate for the remaining individuals and becomes self-enforcing). Preliminary analyses suggest that dampening or amplification of net dispersal flows by synchronous and asynchronous oscillations, respectively, creates different feedback loops based on these different drivers of emigration, but more detailed analyses are required to understand how these contrasting states stabilize themselves.

Effect of periodic environmental disturbances

Periodic environmental disturbances have a stronger effect on population variability on all spatial scales than local interactions or mean patch isolation. We infer this from the observation that both weak and strong trophic cascades, which behave very differently in static landscapes, exhibit almost identical variability patterns in dynamic landscapes, with elevated levels of \(\alpha\)- and \(\gamma\)-variability and low \(\beta\)-variability. Further, all three variability measures are almost constant over a wide range from low-to-medium mean patch isolation. Only at high mean patch isolation, where the patch networks begin to decompose into several isolated components anyway, the effect of the periodic disruption of the patch networks by the blinking patches dwindles and the variability measures become more similar to their values in static landscapes again. Both the increase in \(\alpha\)-variability and the synchronization of the patches, due to the periodic environmental disturbances, are of course not unexpected. The blinking of the patches increases \(\alpha\)-variability by causing low-frequency biomass oscillations through the extinction and recolonization process and by decreasing the mean biomass densities on these patches. Similarly, environmental fluctuations have long been known to be able to synchronize ecological dynamics in coupled habitats (Moran 1953). More surprising is, however, the overruling strength of the effect of periodic environmental disturbances, considering that a blinking cycle (period length \(\lambda =6000\)) is about 150 times slower than the period length of the population cycles in the weak trophic cascade.

Our approach of modeling periodic environmental disturbances as dynamic landscapes, where some patches become periodically uninhabitable, is inspired by the natural example of kettle holes that have a species-rich community during the colder and wetter seasons, but can run dry during the summer (Kalettka and Rudat 2006). Such periodic (in the example: seasonal) environmental disturbances are a common feature of ecological systems, since in most environments seasonally fluctuating climatic drivers exist (Fretwell 1972). Together with the above discussed surprisingly strong effect of even very rarely occurring disturbances, this may explain why empirically observed effects of patch isolation are often small and inconclusive (Fahrig 2003). Environmental disturbances (especially seasonal ones) of course do not always lead to the abrupt extinction of entire local communities, but could, for example, simply modify resource availability or mortality rates. An interesting avenue for future research might therefore be to explore whether such less drastic disturbances also have the potential to overrule the effects of local interactions and landscape configuration. Furthermore, resting stages can play a critical role in the recolonization of periodically uninhabitable patches (Wade 1990). Accounting for them in the model might decrease synchronicity, as they allow for an independent restart of the local communities.

Relevance and effects of dispersal assumptions

Details of the way species dispersal is implemented within a model can have major implications for the arising population dynamics. In nature, a multitude of causes affects an individual’s decision to leave its home patch (Bowler and Benton 2005), among them being, for example, intraspecific competition (Herzig 1995), quality of food resources (Kuussaari et al. 1996) or top-down pressure through parasitism or predation (Sloggett and Weisser 2002). In our model, we use the net growth rate of a species in a given patch to determine its emigration rate. Since the net growth rate depends on both food availability and predation pressure, the model captures multiple of the above-mentioned causes of dispersal. However, we assume that individuals have only knowledge about the growth conditions in the patch they are currently in and not about the conditions in potential target patches. The dispersal rate between any two patches thus only depends on the local conditions in the source patch and on the spatial arrangement of the patches. Using a consumer-resource model with two patches, Abrams and Ruokolainen (2011) showed that when the dispersal rate depends on the difference of the growth rates between source and target patch, asynchronous (antiphase) cycles frequently occur, which promotes stability. With our approach, we only find asynchronous dynamics in static landscapes, but even then synchronous metacommunity dynamics frequently occur.

Conclusions

We conclude that due to indirect effects of local ecological interactions, dispersal is not necessarily a “double-edged sword” (Hudson and Cattadori 1999) (dubbed so because too much of it can synchronize metacommunity dynamics and increase the risk of correlated extinctions), but also that a portfolio effect due to asynchronous oscillations may not always result in reduced variability at the metacommunity level. Furthermore, in each unique landscape, comprising a multitude of abiotic factors, the impact of a periodic environmental disturbance has the potential to outweigh local interactions present in a community. The extent of the effect of mean patch isolation on the variability of population dynamics in a metacommunity thus may strongly depend on local environmental conditions which are relevant for reliable predictions. Whether this is also true for other aspects of fragmentation or habitat loss is an intriguing question for future investigations. Finally, the non-monotonous stability response curve of the strong trophic cascade shows that the effect of mean patch isolation on metacommunity dynamics may not be trivial and that there might be transitions where patch isolation might switch from having a positive to having a negative effect.