Strange invaders increase disturbance and promote generalists in an evolving food web
Model: network structureCommunities are simulated using a modified version of the evolutionary food web models developed in Allhoff et al. (2015) and Allhoff & Drossel (2016), which build on previous models25,26 to show that biodiversity can be maintained in multitrophic networks despite ongoing species turnover when feeding traits are allowed to evolve independent of body mass. The model includes consumptive and competitive interactions, where interaction strengths are determined by the traits of consumer species and their resources. All species possess three traits, a body mass or size ((m)) (used interchangeably), which places them on a body size trait axis, a feeding center ((f)) and feeding range ((s)), which determine the shape and placement of their feeding curve along the axis (Fig. 1a). While the (s) parameter specifically represents one standard deviation of a species’ feeding curve, we refer to (s) throughout as simply the feeding range. The feeding curve represents the hypothetical, fundamental feeding niche of species and shows the potential strength of a consumer’s attack rate for a given resource located along the body size trait axis. Because interactions are determined through these Gaussian curves, our networks are technically fully connected. However, when resources are far from consumer’s feeding centers, interaction strengths become asymptotically small, having a negligible effect on dynamics. Additionally, a basal resource drives energy flow in the food web (Fig. 1a). A summary of all model parameters and variables is provided in Table 2.Model: population dynamicsDynamics are governed by a bioenergetics consumer-resource model, where parameters are scaled to the body mass of species, following previous developments in Yodzis & Innes (1992) and Brose et al. (2006). The rate of change of consumer biomass (({B}_{i})) is given by:$$frac{{dB_{i} }}{dt} = mathop sum limits_{j = resources} e_{j} g_{ij} B_{i} B_{j} – mathop sum limits_{j = consumers} g_{ji} B_{i} B_{j} – mathop sum limits_{j = competitors} c_{ij} B_{i} B_{j} – x_{i} B_{i}$$
(1)
where ({e}_{j}) represents the efficiency of biomass conversion of resource (j) by consumers, ({g}_{ij}) is the mass-specific consumption rate of resource (i) by consumer (j), ({c}_{ij}) is the interference competition between consumer (i) and (j), and ({x}_{i}) is the mass-specific biomass loss from respiration and mortality for consumer (i). The rate of change in basal resource biomass (({B}_{0})) is described by:$$frac{{dB_{0} }}{dt} = n_{0} – mathop sum limits_{j = consumers} g_{j0} B_{j} B_{0} – lB_{0}$$
(2)
where ({n}_{0}) represents the constant influx of resource biomass and (l) the outflow rate. The time scale of the whole system is therefore defined by setting the constant resource influx rate ({n}_{0}=1), meaning that all other rates in the system, and consequently also consumer lifespans, must be interpreted in relation to ({n}_{0}). The basal resource is given a constant body mass trait value of ({m}_{0}=1) which does not evolve. The mass-specific consumption rate is given by:$${g}_{ij} = frac{1}{{m}_{i}} frac{{a}_{ij}}{1+{sum }_{k=prey}{h}_{i}{a}_{ik}{B}_{k}}$$
(3)
where,$${a}_{ij}= {m}_{i}^{0.75}cdot {N}_{ij}={m}_{i}^{0.75}cdot frac{1}{{s}_{i}sqrt{2pi }}cdot mathrm{exp}left[-frac{{left({log}_{10}left({f}_{i}right)-{log}_{10}({m}_{j})right)}^{2}}{2{s}_{i}^{ 2}}right]$$
(4)
describes the mass-specific attack rate of consumer (i) on resource (j), given the feeding kernel (({N}_{ij})) of consumer (i). Gaussian feeding kernels are calculated from consumer (i)’s feeding range (({s}_{i})), feeding center (({f}_{i})), and resource j’s body mass (({m}_{i})), such that resources which occur close to consumer feeding center on the body size trait axis result in the highest attack rates (Fig. 1a). The mass-specific handling time for consumers is given by ({h}_{i}=0.4cdot {m}_{i}^{-0.25}). Interference competition between consumer (i) and (j) is described by:$${c}_{ij}= {c}_{0}cdot frac{{I}_{ij}}{{I}_{ii}} text{ for }ine j$$
(5)
where,$${I}_{ij}= int {N}_{ik}cdot {N}_{jk}dleft({log}_{10}{(m}_{k})right)$$
(6)
describes the overlap in resources (k) between two competing consumers (i) and (j), such that consumers with similar feeding traits will have greater overlap between their feeding kernels resulting in higher competition coefficients.Model: community assembly & network evolutionCommunity assembly of food webs occurs through a combination of ecological and evolutionary dynamics (Fig. 1b). All ecological dynamics are described by the consumer-resource model above, where species with viable biomass densities persist in communities and species whose biomass falls below a fixed extinction threshold ((varepsilon = {10}^{-8})) are removed from the network. New species are introduced probabilistically into the network at fixed intervals through either mutation events ((p)) or as invaders ((1-p)), where (p) can be manipulated to increase the frequency of either mutation or invasion events. The traits of new mutant species are drawn probabilistically from a Gaussian distribution set around the traits of a selected extant parent species in the network. Invader species traits are generated in a similar fashion but using a Gaussian distribution with a greater standard deviation. The standard deviation of this trait range is set with the invader strangeness parameter (z), which can be manipulated to increase the range of potential traits for invader species. Thus, a larger (z) value increases the probability that new invader species will appear “strange” compared to other species already in the community. For mutant species, (z) is always set to 0.1.Parents of mutants are chosen probabilistically, where species with greater individual density (species biomass/body mass) are more likely to generate new mutant species. The parents of invader species are chosen randomly, with equal probability given to all extant species in the community. Both mutants and invaders are introduced into the system at the extinction threshold biomass ((varepsilon ={10}^{-8})). For mutants the initial biomass is removed from the biomass of the parent species’ populations, while for invaders this biomass is added into the system without affecting the parent species’ biomass pool; however, this difference did not significantly impact our results.Communities are initialized with a single ancestor species (starting biomass (varepsilon ={10}^{-8})) and the basal resource (starting biomass (=frac{{n}_{0}}{l}=2.0)) (Fig. 1b). The ancestor species is given a body mass of (m=100), feeding center of (f=1), and feeding range of (s=0.4). Upon initialization, the system is a run with only the ancestor species consuming the basal resource until a new species is introduced at 100 time steps. Thereafter, new species are introduced every 100 time steps, with ecological dynamics occurring between each species introduction. Additionally, species biomass is assessed at each 100 time step interval and non-viable species populations that fall below the biomass extinction threshold are removed. This process is repeated cyclically over the course of simulations (Fig. 1b), with many new species being generated and many removed due to extinction. The persistence of individual species is thus determined by their individual traits and overall resource availability given the composition of the rest of the community. With this dynamical approach to simulating evolving food webs, similar models have been shown to generate viable communities with both multi-trophic diversity and constant species turnover27,28, making this framework useful for testing the evolutionary impacts of species invasion and disturbance on community composition.Simulation experimentsSimulations were conducted in C, where numerical integration of differential equations was performed using the Runge–Kutta–Fehlberg algorithm from the GNU Scientific library29. Simulations were run for 25 million time steps, with 250,000 novel species introductions (mutants or invaders) for each simulation. To test if invasion would increase disturbance and variability in communities and drive the evolution of more generalized species, we conducted simulations where invaders were introduced with an increased probability of having trait values that were divergent from parent species. We controlled this by manipulating the invader strangeness parameter ((z)) across a range from (z=0.1) (invader and mutant trait values are equivalent) to (z=5.0). Invasion frequency ((p)) was fixed at 0.2 for all simulations, making mutation events more likely to occur than invasion.We hypothesized that introducing invaders with traits that are very different from parent species and from the community should result in greater disturbance in food webs because these species would be more likely to occupy novel niche space along the body size trait axis, which could result in the overexploitation of resources either through superior feeding strategies or by allowing invaders to avoid consumption by other consumers. Together, this should increase the probability of disrupted consumer-resource dynamics and secondary species extinction occurring with the introduction of strange invaders, both resulting in increased variability of biomass in the community. As a result, this increased variation should favor the survival of more generalist species in the community if they can buffer variability by consuming a greater range of resources.This is tested against the assumption that specialist consumers are more efficient than generalist consumers (generalist trade-off hypothesis2,12), which is built into our model given the formulation of the attack rate parameter ((a)), where specialist species achieve higher optimal peaks in attack rates, given their smaller feeding ranges ((s)). Thus, under conditions of low variability, our model results in communities being composed of mostly very specialized species, with narrow feeding ranges. To counter this trend toward extreme specialization, we set a floor for minimum feeding range values for all species of (s=0.3). Given these tendencies, we expected the persistence (lifespan) of more specialist species to be greatest under conditions of low variability (low invader strangeness) and that the relative persistence of more generalist species compared to specialists should increase with disturbance due to increasingly strange invaders.To test the robustness of these predictions, we replicated the (z) parameter sweep 100 times using random initial seed sets, resulting in 5000 simulations total, which collectively generated over 1.25e+09 unique species across all simulations. Data from these simulations was extracted at three different time intervals. We assessed species traits and lifespan data for all species generated in simulations at every 100 time steps, excluding data from the first 50,000 time steps to avoid including transient dynamics. Community level data, including community biomass and basal resource biomass were extracted at every 50,000 time steps (excluding time 0 from analysis). Species turnover data was extracted at every 10,000 time steps. In the infrequent event that simulations did not complete (community level extinction or crashed runs) we reran simulations with different random seed sets but identical parameter values.Data & statistical analysisDo resource and community variation increase with invader strangeness?To assess whether the addition of increasingly strange invaders into food web communities resulted in increased variation we analyzed several metrics of community and resource variability. We calculated the standard deviation (SD) of the basal resource biomass across time for each simulation and pooled these data for all simulation replicates across the invader strangeness parameter sweep. To assess variation at the community level, we used a similar approach to calculate variability in community biomass. For this metric, we summed the population biomasses of all species in the community for each given time interval output (excluding species introduced at that time step) and calculated the SD of these values across time for each individual simulation.Finally, to further assess community variability and to determine if increasing invader strangeness drives increased extinction in communities, we calculated species turnover for each time output. Species turnover was measured as the percentage change in the composition of species in communities between each time output (10,000 time steps). We then calculated mean species turnover over time for each simulation replicate and pooled all data together. To account for the non-linearity observed in our variation data (see “Results”) we conducted generalized additive models (GAM) to determine if increasing invader strangeness resulted in a significant increase in variability. GAMs were fit using a gamma error distribution with a log link function to account for continuous data constrained to positive values.Does the degree of generalism in communities increase with invader strangeness?To determine if the degree of generalism and the proportion of generalist species in food webs increased as invader strangeness increased, we calculated the mean and median feeding range ((s)) (Table 1) of species which occurred in communities for each simulation. We included all species that were generated and that survived for at least 100 time-steps in simulations, to remove the many non-viable species which immediately go extinct. Additionally, we included only mutant species for this metric to avoid the influence of the traits of invaders species, which we directly manipulated through the invader strangeness parameter. We reasoned this would provide a more independent metric of feeding range trends in communities. Mean and median feeding range were calculated for all simulation replicates and the impact of invader strangeness was assessed with GAMs (gamma error distribution with a log link function) to account for non-linear data (see “Results”).Additionally, we calculated a measure of the realized feeding range of consumers (distinct from the fixed fundamental feeding range ((s)) (Table 1)) to determine if more species were functioning as feeding generalists in communities. For this metric, we calculated the attack rate of each consumer on all other species in the community (including the basal resource and the focal consumer) for each time output (every 50,000 time steps from our community data, excluding species introduced at that time step). We then calculated the proportion of the attack rate on each species compared to the focal species’ maximum possible attack rate (an ideal prey at the exact center of the consumer’s feeding kernel). We then excluded all values below a threshold of 0.1 and from this calculated the proportion of species consumed out of the total number of species in the community. This metric correlated positively with the fundamental feeding range ((s)) of consumer species (Supplementary Fig. S3) and we refer to it throughout as the realized feeding range (Table 1) of consumer species. For our statistical analysis, we calculated the mean realized feeding range of species per simulation across invader strangeness ((z)) and ran a GLM with a quasibinomial error distribution and logit link function to account for proportional data.Does the persistence of generalist species increase with invader strangeness?To determine the persistence of species in our simulations we assessed the lifespan of individual species in simulated communities across time. For a given species, lifespan was measured as the number of time steps it persisted in a simulation after its initial introduction. We used this data to determine the relationship between species persistence and feeding range traits in two ways. First, we assessed the lifespan of all species in individual simulations continuously given the feeding range trait values across species. From this, a regression coefficient was calculated from the log10-scaled data, using a GLM with a gamma error distribution (log link), to determine the trend or “lifespan slope” for each simulation under different levels of invader strangeness (Fig. 4b). These lifespan slope values were then assessed for all simulation replicates across the full range of the invader strangeness parameter. Because more specialized species have higher maximal attack rates and are typically more efficient in our model, we expected that the lifespans of specialist species would be longer than more generalized species and that lifespan slopes should be negative under conditions of low variation. Given this, we expected to observe a positive trend in lifespan slope values across the invader strangeness parameter sweep if disturbance was increased in simulations as (z) became higher. We tested for this positive trend in the lifespan slope data by conducting a GAM (Gaussian distribution and the identity link function) to manage the observed non-linear trend in our data (see “Results”).For the second approach, we aimed to determine the relative persistence of species by binning “generalist” and “specialist” species based on feeding range traits and comparing species lifespans between these groups. For this analysis we split species into bins, where specialist species included all species with feeding range (sle) 0.32 and generalists as all species with feeding range values (sge) 0.39 (species with intermediate feeding range values were excluded from the analysis). We performed a robustness check of bin cutoffs but found no qualitative or statically significant differences in our results for a range of bin cutoff values. To assess how the relative persistence of generalists compared to specialists was influenced by invader strangeness, we then calculated the mean life span of all species falling into either of these categories per simulation and determined how these values were influenced by (z) for all simulation replicates. To assess whether mean lifespan was different between each of these groups across the invader strangeness sweep, we conducted a GLM with species type (generalist or specialist) and invader strangeness ((z)) as fixed effect terms and tested for the statistical significance of their interaction on mean species lifespan. The GLM was run using a Poisson distribution to account for discrete lifespan count data with a log link function. All GLMs and GAMs were performed in R using the “glm” and “mgcv” functions30, respectively, and all non-linear parameters in GAMs were fit using generalized cross validation (GCV). More
