Landscape heterogeneity buffers biodiversity of simulated meta-food-webs under global change through rescue and drainage effects

Model

We model a tri-trophic food chain of one plant, one herbivore and one predator population on one or two habitat patches and complex meta-food-webs consisting of 10 plants and 30 animals in different landscapes containing 50 patches. The feeding dynamics are constant overall patches and are determined by the allometric food-web model by Schneider et al. 2016³³. We integrate dispersal as species-specific biomass flux between habitat patches according to Ryser et al. 2019³⁴. With the use of a dynamic bioenergetic model we formulate feeding and dispersal dynamics in terms of ordinary differential equations. The rate of change in biomass densities of a species are the sum of its biomass loss by metabolism, being preyed upon and emigration and its biomass gain by feeding and immigration. For detailed equations and for model parameters see section Equations and parameters and the supplement (Supplementary Table 1).

Local food-web dynamics

Following the allometric food-web model by Schneider et al. 2016³³ each species is fully characterised by its average adult body mass. For the complex food-web log₁₀ body masses were randomly drawn from a uniform distribution from 0 to 3 for plants and from 2 to 6 for animals. For the food chain the plant body mass was set to 10², the herbivore body mass to 10⁴ and the predator body mass to 10⁶. We set mass ratios of the herbivore to the plant and the predator to the herbivore to the optimum of 100, thus the respective resource being a one-hundredth of its consumer’s body mass. This simplifies feeding efficiency rates (see section Equations and parameters; L_i,j, Eq. (5)) to 1 in the case of a food chain. Trophic dynamical parameters, such as metabolic rates and feeding rates, scale with body masses of model species. Also, we assume a type-II functional response for the food chain and a slight nonlinearity of the functional response in the food web as this stabilises persistence in more complex systems. Compared to Ryser et al 2019, capture rates were reduced to 5% to achieve viable food chains and food webs to increase the stability in the absence of interference competition.

Nutrient model

We have an underlying nutrient model with one nutrient that is driving the nutrient uptake and therefore the growth rate of the plant population^11,33. The nutrient model consists of one nutrient, a nutrient turnover rate of 0.25 and a nutrient supply concentration. The nutrient supply concentration was varied to get eutrophic and oligotrophic patches (see Setup).

Spatial dynamics

We model dispersal between local communities as a dynamic process of emigration and immigration, assuming dispersal to occur at the same timescale as the local population dynamics³⁵. Thus, biomass flows change dynamically between local populations and the dispersal dynamics directly influence local population dynamics and vice versa²⁵.

Dispersal rates of animals are modelled with an adaptive emigration rate depending on the net growth rate on the given patch. Dispersal ranges depend on the body masses of our model species with larger species having a higher dispersal range. We model a hostile matrix between habitat patches that does not allow feeding interactions to occur during dispersal. Depending on the scenario, we define a landscape with one, two or 50 patches. In cases with two or 50 patches, their locations are spatially explicit and were chosen in a way that the distances between reflect the dispersal loss of the predator across the matrix hostility gradient.

Emigration and immigration

Based on empirical observations³⁶ and previous theoretical frameworks^13,22,37, we assume that the maximum dispersal distance of animal species increases with their body mass. For simplicity, we do not let the plants disperse, as they do not move themselves and the dispersal of plant propagules strongly depends on their dispersal strategy. We model emigration rates as a function of each species’ per capita net growth rate, which is summarising local conditions such as resource availability, predation pressure, and inter- and intra-specific competition²⁵ (but see Sensitivity Analyses for dispersal models with constant dispersal or non-body-mass-scaled dispersal ranges). Dispersal losses scale linearly with the distance between two patches and are 100% in scenarios with only one patch or when the distance between the two patches surpasses the dispersal range of an animal. Even though we model dispersal losses according to dispersal distances, this loss term could also represent any other sort of dispersal loss. For numerical reasons, we did not allow dispersal flows smaller than 10⁻¹⁰.

Numerical simulations

We initialised each local population with a biomass density randomly sampled from a uniform probability density within the interval (0,10). Starting from these random initial conditions, we numerically simulated food web and dispersal dynamics over 100,000 time steps by integrating the system of differential equations implemented in C++ using procedures of the SUNDIALS CVODE solver version 2.7.0 (backward differentiation formula with absolute and relative error tolerances of 10⁻¹⁰) and the time series of biomass densities were saved for last 10,000 time steps. For numerical reasons, a local population was considered extinct and was set to 0 once its biomass density dropped below 10⁻²⁰. Based on the empirically derived metabolic rates, these 100,000 time steps correspond to ~11 years. Our model does, however, not account for time spent for organisms’ other non-trophic activities such as sleeping or mating. Thus, the time scales of the simulation should only be compared with caution to natural time scales of population dynamics. Transient dynamics usually equilibrate within the first few thousand time steps.

Equations and parameters

Our model formulates the change of biomass densities over time in ordinary differential equations. Given the empirical origin of metabolic rates used in our model, one time step corresponds to an hour and body masses are in mg, areas of patches are not defined. The feeding links (i.e. who eats whom) are constant overall patches and are as well as the feeding dynamics determined by the allometric food-web model by Schneider et al. 2016³³. We integrate dispersal as species-specific biomass flow between habitat patches. Using ordinary differential equations to describe the feeding and dispersal dynamics, the rate of change in biomass density B_i,z of species i on patch z is given by

with the first three terms describing local trophic dynamics and the last two terms describing emigration, E_i,z (Eq. 9), and immigration, I_i,z (Eq. 11). For simplicity, we do not let plants disperse. Trophic dynamics are driven by following three processes. First, predation or herbivory on species j with assimilation efficiency e (e_j = 0.545, if j is a plant, typical for herbivory; e_j = 0.906 if j is an animal, typical for carnivory³⁸) and the functional response F_ij,z (Eq. 3) for animals, and a nutrient dependent growth (Eq. 7) for plants. Second, losses due to predation or herbivory, respectively. Third, losses by metabolic demands with x_i = x_Am_i^−0.305 with scaling constant x_A = 0.141 (tenfold laboratory metabolic rate³⁹ at a temperature of 20° Celsius to represent field metabolic rates) for animals and x_i = x_Pm_i^−0.25 with x_P = 0.138 for plants. We used a dynamic nutrient model (Eq. 8) as the energetic basis of our food web. Each species i is fully characterised by its average adult body mass m_i. Body masses determine the interaction strengths of feeding links as well as the metabolic demands of species. Data from empirical feeding interactions are used to parametrise the functions that characterise the optimal prey body mass and the location and width of the feeding niche of a predator³³. From each m_i a unimodal attack kernel, called feeding efficiency L_ij is constructed which determines the probability of consumer species i to attack and capture an encountered resource species j. We model L_ij as an asymmetrical hump-shaped Ricker’s function (Eq. 5) that is maximised for an energetically optimal resource body mass (optimal consumer-resource body mass ratio R_opt = 100) and has a width of γ. The maximum of the feeding efficiency L_ij equals 1. Supplementary table 1 is an overview of the standard parameter set for the equations. See also Schneider et al. 2016³³ for further information regarding the allometric food-web model.

Functional response

Per unit biomass feeding rate of consumer i as function of the biomass density of the resource R_j, with b_i,j, resource-specific capture coefficient (Eq. 4); h_i,j, resource-specific handling time (Eq. 6); ω_i = 1/(number of resource species of i), an inefficiency parameter for generalists assuming that generalist are less adapted in for example search patterns or hunting strategies to a specific prey species; and q, the Hill coefficient for nonlinearities in density dependency (if q = 0 it is a Type-II functional response, if q = 1 it is a Type-III functional response).

Capture coefficient

Resource-specific capture coefficient of consumer species i on resource species j scaling the feeding kernel L_ij by a power function of consumer and resource body mass, assuming that the encounter rate between consumer and resource scales with their respective movement speed. This body mass scaling of encounter rates is assumed to occur before the attempt of a predator to capture its prey is made. We differentiate between carnivorous and herbivorous interactions with each comprising a constant scaling factor for their capture coefficients a_k with k ∈ 0, 1 (a₀ = 15 for carnivorous species and a₁ = 3500 for herbivorous species). For plant resources, ({m}_{j}^{{beta }_{j}}) was replaced with the constant value of 1 (as plants do not move).

Feeding efficiency

The probability of consumer i to attack and capture an encountered resource j (which can be either plant or animal), described by an asymmetrical hump-shaped curve (Ricker’s function), centered around an optimal consumer-resource body mass ratio R_opt = 100³³ and with γ that that affects the width of the hump. An increase in γ results in a decrease in the width.

Handling time

The time consumer i needs to kill, ingest, and digest resource species j, with scaling constant h₀ = 0.4 and allometric exponents η_i = −0.48 and η_j = −0.66⁴⁰.

Growth factor for plants

Species-specific growth factor of plants determined dynamically by the nutrient; with K_i, half-saturation densities determining the nutrient uptake efficiency assigned randomly for each plant species i and (uniform distribution within (0.1, 0.2)).

Nutrient dynamics

Rate of change of nutrient concentration N of nutrient on patch z, with global turnover rate D = 0.25, determining the rate at which nutrients are refreshed and the nutrient supply concentration S.

Generating landscapes

We generated different fragmented landscapes, represented by random geometric graphs, by randomly drawing the locations of Z patches from a uniform distribution between 0 and 1 for x- and y-coordinates, respectively.

Dispersal

We model dispersal between local communities as a dynamic process of emigration and immigration, assuming dispersal to occur at the same timescale as the local population dynamics. Thus, biomass flows dynamically between local populations and the dispersal dynamics directly influence local population dynamics and vice versa. We model a hostile matrix between habitat patches that does not allow for feeding interactions to occur during dispersal. The total rate of emigration of animal species i from patch z is

with d_i,z as the corresponding per capita dispersal rate. We model d_i,z as

with a, the maximum dispersal rate, b = 10, a parameter determining the shape of the dispersal rate, x_i, the inflection point determined by the metabolic demands per unit biomass of species i, and υ_i,z, the net growth rate of species i on patch z. The net growth rate consists of the biomass gain by feeding, the biomass loss by being fed upon and the metabolic loss (({v}_{i,z}=frac{{B}_{i,z}mathop{sum}limits_{j}{e}_{j}{F}_{{ij},z}-mathop{sum}limits_{j}{{B}_{j,z}F}_{{ji},z}-{x}_{i}{B}_{i,z}}{{B}_{i,z}})). We chose to model d_i,z as a function of each species’ net growth rate to account for emigration triggers, such as resource availability, predation pressure, and inter- and intra-specific competition. If for example an animal species’ net growth is positive, there is no need for dispersal and emigration will be low. However, if the local environmental conditions deteriorate, the growing incentives to search for a better habitat increase the fraction of individuals emigrating.

Immigration

The rate of immigration of biomass density of species i into patch z follows

where N_z and N_n are the sets of all patches within the dispersal range of species i on patches z and n, respectively. In this equation, E_i,n is the emigration rate of species i from patch n, ({max }(1-{delta }_{i,{nz}},0)) is the fraction of successfully dispersing biomass, i.e. the fraction of biomass not lost to the matrix, and δ_i,nz is the distance between patches n and z relative to species i’s maximum dispersal distance δ_i (see below paragraph Maximum dispersal distance). The term (frac{{max }(1-{delta }_{i,{nz}},0)}{mathop{sum}limits_{m,epsilon, {N}_{n}}{max }(1-{delta }_{i,nz},0)})determines the fraction of biomass of species i emigrating from source patch n towards target patch z. This fraction depends on the relative distance between the patches, δ_i,nz, and the relative distances to all other potential target patches m of species i on the source patch n, δ_i,nm. Thus, the flow of biomass is greatest between patches with small distances to account for the logic that the first patch dispersing organism come across is closer. In other words, the further a destination is, the more likely it is to come across another patch before.

For numerical reasons, we did not allow for dispersal flows with I_i,z < 10⁻¹⁰. In this case, we immediately set I_i,z to 0. We assume that the maximum dispersal distance δ_i of animal species increases with their body mass. For animal species, the body mass m_i determines how far they can travel through the matrix. Thus, animal species at high trophic positions can disperse further than smaller animals at lower trophic levels. Each animal species perceives its own dispersal network dependent on its species-specific maximum dispersal distance

where the exponent ε = 0.05 determines the slope of the body mass scaling of δ_i and intercept δ₀ = 0.1256. This intercept had been chosen because an organism with a body mass of 10¹² would have a maximum dispersal range of 0.5. We chose a positive value for ε to account for a higher mobility of animals with larger body masses.

Setup

To answer our questions, we model the following scenarios:

Nutrient enrichment (Fig. 2a): Simulations across a gradient of nutrient supply concentrations (0, 10) on one patch without emigration and therefore also no dispersal loss.

Drainage effect (Fig. 2b): Simulations across a gradient of maximal emigration rates (0, 0.15) on one eutrophic patch with a nutrient supply concentration of 10.

Hostility effect with two patches: Simulations across a gradient of dispersal losses (0, 1) on two eutrophic patches with a nutrient supply concentration of 20 on each and a maximal dispersal rate of 0.05.

Heterogeneity effect with two patches (Fig. 3, along the x-axis): Simulations across a gradient of nutrient supply concentrations (0, 20) on one of two patches with the other patch being a eutrophic patch with a nutrient supply concentration of 20, a maximal emigration rate of 0.05 and no dispersal loss.

Interaction of hostility effect and heterogeneity effect (Fig. 3): For each level of heterogeneity (difference in nutrient supply between the two patches) we simulated the whole gradient of the hostility effect (dispersal loss of the predator from 0 to 1).

Heterogeneity effect on complex food webs in complex landscapes (Fig. 4): For a complex meta-food-web, we generated five random geometric graphs consisting of 50 patches. Each patch was initialised with a complex food web consisting of 10 plant and 30 animal species. For each random geometric graph, we simulated 15 homogeneous landscapes, where all patches have the same nutrient supply concentration with simulations across a gradient of nutrient supply concentrations ranging from 10^−0.8 (oligotrophic) to 10² (eutrophic) in steps of 0.2 in the exponent, and 25 heterogeneous landscapes for each of the three heterogenous scenarios (heterogeneous oligo-, meso-, and eutrophic), where the nutrient supply concentration for each patch is assigned randomly from the same gradient as in the homogeneous scenario. For the heterogeneous mesotrophic landscapes, the gradient from 10^−0.8 to 10² (in steps of 0.2 in the exponent) was sampled from a uniform distribution, while for the heterogenious oligotrophic and eutrophic landscapes, 42 nutrient supply values were sampled from the lower (10^−0.8 to 10⁰) or higher third (10^1.2 to 10²) of the gradient, repectively, 5 from the intermediate third (10^0.2 to 10¹) and 3 of the remaining third of the gradient.

Sensitivity analyses

Additional results from simulations with a non-adaptive dispersal model and non-body mass-scaled dispersal ranges of organisms are presented in the Supplementary Notes. Results are, however, qualitatively the same.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com