Modeling foundation species in food webs

. Foundation species are basal species that play an


INTRODUCTION
Foundation species (sensu Dayton 1972) are basal species that structure ecological communities by creating physical structure and modulating ecosystem processes (Ellison et al. 2005).Recent declines (e.g., Tsuga canadensis) and extirpations (e.g., Castanea dentata) of foundation species in terrestrial ecosystems have called attention to the need for new methods for identifying and quantifying the role of foundation species in ecological communities (reviewed by Ellison et al. 2005;2010, Van der Putten 2012).Numerous field studies have shown that foundation species can alter trajectories of the assembly of ecological communities (e.g., Gibson et al. 2012, Schoeb et al. 2012, Butterfield et al. 2013, Martin and Charles 2013, Orwig et al. 2013).However, general models of how foundation species affect ecological systems are scarce and generally qualitative (Ellison and Baiser, in press).
Foundation species can interact trophically within a community, but they exert their influence primarily through non-trophic effects (Ellison and Baiser, in press).Some examples of non-trophic actions of foundation species include; altering local climates and microclimates (e.g., Schoeb et al. 2012, Butterfield et al. 2013); changing soil temperature, moisture, and acidity (e.g., Prevey et al. 2010, Lustenhouwer et al. 2012, Martin and Charles 2013); providing refuge for prey species and perches for predators (e.g., Yakovis et al. 2008, Tovar-Sanchez et al. 2013) ; and stabilizing stream banks and shorelines against erosion (reviewed by Ellison et al. 2005).
Because foundation species exert system-wide effects on biodiversity and ecosystem functioning primarily through these (and other) non-trophic interactions, it has proven difficult to link effects of foundation species into theories of the structure and function of food webs.Food-web theory aims to elucidate the persistence of the types of complex, species-rich webs that we see in nature (e.g., May 1972, Allesina andTang 2012).Measures of network properties, such as connectance, compartmentalization, and species richness, as well as the strength of species interactions, all can influence the stability and persistence of food webs (e.g., May 1972, Dunne et al. 2002, Gravel et al. 2011, Stouffer and Bascompte 2011).Adding non-trophic interactions, such as those exhibited by foundation species or mutualists in general, provides an additional step towards understanding persistence and stability of ecological networks (Thébault and Fontaine 2010, Allesina and Tang 2012, Kéfi et al. 2012) Here, we adapt non-linear, bioenergetic predator-prey models to explore non-trophic roles of foundation species in food webs.To make explicit linkages between trophic and nontrophic interactions, we model the metabolic rate of individual "species" as a function of foundation species biomass.Metabolic rate is good proxy for a wide variety of positive nontrophic species interactions (sensu Kéfi et al. 2012), because "stressful conditions" may be reduced when foundation species ameliorate temperature extremes, provide associated species with habitat resources or shelters, or enhance their growth rate (Schiel 2006, Shelton 2010, Gedan et al. 2011, Angelini and Silliman 2012, Dijkstra et al. 2012, Noumi et al. 2012, Butterfield et al. 2013).
We developed four different foundation species models to explore non-trophic effects of foundation species in food webs.In each, the foundation species influences target species at different trophic positions in the food web: 1) a basal model, in which the foundation species reduces the metabolic rates of only other, albeit non-foundation, basal species, 2) a consumer model, in which the foundation species reduces the metabolic rates of only consumers, 3) a total model, in which the foundation species reduces the metabolic rates of all species, and 4) a control model, in which the foundation species is only consumed and has no effect on the metabolic rates of any associated species.We examined the outcomes of each of these models for six metabolic rate "treatments" in which the foundation species alters the metabolic rates of associated species by one-tenth to ten times their allometric baseline metabolic rates.For each model simulation, we looked at how foundation species influence different measures of foodweb structure during community assembly and the subsequent change of food-web structure when the foundation species was removed.

METHODS
We modeled dynamic ecological networks using a four-step process (Brose et al. 2006, Berlow et al. 2009, Kéfi et al. 2012): 1) model initial network structure; 2) calculate body mass for each species based on trophic level; 3) simulate population dynamics using an allometric predator-prey model; and 4) add non-trophic interactions into the allometric predator-prey model.

Network structure
We used the niche model of Williams and Martinez (2000) to designate trophic links in our model food webs.The niche model is an algorithm with two parameter inputs: species richness (S) and connectance (C = L/S 2 , where L = the number of trophic links).Each species in the web has a niche value uniformly drawn from [0,1] and a niche range that is placed on a onedimensional axis.Any one species whose niche value falls within the niche range of another is defined to be the latter's prey (for specific details on the niche model see Williams and Martinez 2000).The niche model has been shown to reproduce accurately a wide range of food-web network properties for many empirical webs (Williams and Martinez 2000, Dunne et al. 2004, Williams and Martinez 2008).

Body mass
We calculated body mass, M i , for species i as: In eq. ( 1), Z is the predator-prey biomass ratio and T is the average trophic level of species i calculated using the prey-averaged method (Williams and Martinez 2004).We set basal species M to unity and used a predator-prey biomass ratio of Z = 10 2 .We used body mass to allometrically scale biological parameters in the predator-prey model.

Allometric predator-prey model
We simulated food-web population dynamics using an allometric predator-prey model (Yodzis and Innes 1992, Williams and Martinez 2004, Brose et al. 2006).Following Brose et al. (2006): Equation 2a describes change in biomass, B, of primary producer species i, and equation 2b describes changes in B of consumer i.All model variables are listed and defined in Table 1.
For primary producer species i, r i is its mass-specific maximum growth rate; M i is its individual body mass; and G i is its logistic growth rate: G i = 1-(B i /K) and K is the carrying capacity (in our model, K = 1).Both for primary producers and consumers, the mass-specific metabolic rate for species i is x i .For consumers, y i is the maximum consumption rate of species i relative to its metabolic rate; e ji is the assimilation efficiency for species i when consuming species j; and f ij is the fraction of biomass removed from the resource biomass that is actually ingested.The functional response, F ij , describes how consumption rate varies as a function of prey biomass.We used a type II functional response: In eq. ( 3), ω ij is the uniform relative consumption rate of consumer i preying on resource j (i.e., the preference of consumer i for resource j) when the consumer has n total resources (ω ij = 1/n) and B 0 is the half-saturation constant (i.e., resource biomass at which consumer reaches half of its maximum consumption rate).In all of our models, B 0 was set equal to 0.5.
Body size is an important component of both predator-prey interactions (Warren and Lawton 1987, Woodward and Hildrew 2002, Brose et al. 2006) and metabolic functioning of organisms (Brown et al. 2004).As a result, body size is an important factor for energy flow throughout food webs (Woodward et al. 2005).Predator-prey body-size ratios found in empirical food webs have been shown to stabilize dynamics in complex networks (Brose et al. 2006).
Thus, we allometrically scaled the biological parameters r i , x i , and y i in eqns (2a) and (2b) to body size (Brose et al 2006).We modeled the biological rates of production, R, metabolism, X, and maximum consumption rate, Y, using a negative-quarter power-law dependence on body size (Brown et al. 2004): In eqns (4a-4c), subscripts P and C correspond to producers and consumers respectively; a r , a x , and a y are allometric constants; and M is the body mass of an individual (Yodzis and Innes 1992).The time scale of the system is specified by fixing the mass-specific growth rate, r i , to unity.Following this, we normalized the mass-specific metabolic rate, x i , for all species in the model by time scale and in turn, we normalized the maximum consumption rate, y i , by the metabolic rates: We then entered the allometrically scaled parameters for r i , x i , and y i into equations 2a and 2b, yielding an allometrically scaled, dynamic predator-prey model.We set the allometric constants to be y i = 8, e ij = 0.85 for carnivores and e ij = 0.45 for herbivores, a r = 1, and a x = 0.314 (Yodzis and Innes 1992, Brown et al. 2004, Brose et al. 2006).

Foundation species and non-trophic interactions
For each food web, we randomly designated one basal species as a foundation species.
Each foundation species engaged in a non-trophic interaction with a given number of target species in a food web, depending on the model described in the next section.The foundation species alters the metabolic rate (x) of a target species with which it interacts following a general saturating function (after Otto and Day 2007): In eq. ( 6), x fsp is the metabolic rate of the target species in the presence of the foundation species; x a is the metabolic rate of the target species in the absence of the foundation species (i.e., baseline metabolic rate, eqn 5b); B is the biomass density of the foundation species; and B a is the "typical" (i.e., ~average across trial runs) biomass density for the foundation species.The metabolic rate of species i, x i , decreases from x a when B = 0 to an asymptote at x fsp when B is large (we assume that x fsp < x a because the foundation species reduces the metabolic rates of its associated species).

Four foundation species models
We varied the number and position of non-trophic interactions in four different ways (Fig. 1).In the control model, there are no non-trophic interactions (i.e., the species designated as the foundation species has only trophic interactions).In the basal model, the foundation species influences the metabolic rate of all basal species.In the consumer model, the foundation species influences the metabolic rate of all consumers (i.e., non-basal species).Finally, in the total model, the foundation species influences the metabolic rate of all species in the food web.

Simulations and analysis
We created 100 niche-model webs, in all of which we set S = 30 and C = 0.15.We parameterized allometric predator-prey models with an initial biomass (B i ) vector drawn randomly from a uniform distribution: B i ~ Uniform[0.5,1].The initial value of B i was the same for any given food web in all four of the foundation species models.We solved equations 2a and 2b using the standard 4 th order Runge-Kutta method with a time step of 0.001.For each model run, we ran the initial "food-web assembly" simulations for 2,000 time steps.A species was considered extinct and removed from model simulations (i.e., B i = 0) when B i <10 -30 (Brose et al. 2006, Berlow et al. 2009).At the end of this "assembly" period we calculated the number of species present and nine additional measures of food-web structure (Table 2) and then removed food webs with unconnected species or chains from further simulation.We next "removed" the foundation species from the remaining webs and ran the "foundation species removal" simulation for an additional 2,000 time-steps.At t =4,000, we again calculated the number of species present and the nine additional measures of food web structure (Table 2).
Food-web metrics (Table 2) were calculated using Network 3D (Williams 2010).For the food web assembly analysis (i.e., the first 2,000 time steps of each model run), we tested the effect of each model (foundation species effects) using analysis of covariance (ANCOVA).In the ANCOVA, foundation species model was the factor, and log (metabolic rate +1) was the covariate.Because measures of food-web structure are often correlated (Vermaat et al. 2009), we used principle components analysis (prcomp in R version 2.13.1) to reduce the food-web metrics into two orthogonal principle components that were used as response variables in the ANCOVA.
In this analysis, we did not include food webs that collapsed (i.e., had zero species).ANCOVA was implemented using glm in R; a Poisson link function was used when species richness was the response variable, and a Gaussian link function was used for the analysis of food-web metrics (principal axis scores).
For the foundation species removal analyses (i.e., time steps 2,001 -4,000), we calculated standardized change (∆z = z t=2001 -z t=4000 / z t=2001 ) in species richness and food-web metrics (principal axis scores) between the end of food-web assembly (t = 2,001) and the end of the foundation species removal (t = 4,000) because webs had different species richness at the time the foundations species was removed (t = 2,000).As described above, we then used ANCOVA to test the effects of each model.

Exploring the parameter space
An important assumption in our models is that species have higher metabolic rates in the absence of the foundation species.However, it was not clear how to set the baseline metabolic rate, x a , (i.e., how poorly should any particular species perform in the absence of the foundation species) and how much the foundation species should improve [= reduce] the metabolic rate (x fsp ).To explore a range of reasonable possibilities, we ran one set of simulations in which x a was set equal to the allometrically scaled metabolic rate in eqn.(5b) and x fsp was set equal to one of 0.5, 0.2 or, 0.1 of x a (Fig. 2A; referred to henceforth as 0.5×, 0.2×, and 0.1× treatments).In this first set of simulations, species start at the (allometric) baseline and the presence of the foundation species further reduces the metabolic rates of species associated with it.In the second set of simulations, x fsp was set equal to the allometrically scaled metabolic rate in eqn.(5b) and x a was set equal to one of 2, 5, or 10 times x fsp (Fig. 2B; referred to henceforth as 2×, 5×, and 10× treatments).Our metabolic rates encompass the variation observed between basal metabolic rates and maximum metabolic rates in empirical studies (Nagy 1987, Gillooly et al. 2001).
In total, we simulated 100 webs for each combination of the four foundation species models and the six metabolic treatments: 100 × 4 × 6 = 2,400 food-web simulations.Model code is available from the Harvard Forest Data Archive (http://harvardforest.fas.harvard.edu/dataarchive),dataset HF-211.

FOOD-WEB STRUCTURE
The first two principal components of food-web structure (Fig. 4) accounted for 67% of the variation across model food webs (Table 3).Model webs with low PC-1 scores were relatively species-rich with high C, LS, and cluster coefficients, and also had a high fraction of intermediate species and omnivores.Conversely, webs with high PC-1 scores were species-poor with low C and LS; these webs also had long path lengths and large fractions of top, basal, and herbivore species.Webs with high PC-2 scores were species-rich with low C, and had large proportions of top species, low proportions of basal species, and low cluster coefficients.Webs with low PC-2 scores were species-poor with high C and cluster coefficients, and had a large fraction of basal species.

FOOD-WEB STRUCTURE
The first two principal components accounted for 60% of the variation in food-web structure after the removal of the foundation species (Table 3).Model webs with high PC-1 scores lost a greater proportion of species, and showed relatively larger decreases in LS and cluster coefficients (Fig. 6).These structural changes were due primarily to a decrease in the proportion of intermediate and omnivore species and an increase in the proportion of basal species after foundation species removal.Webs with low PC-1 scores lost fewer species and experienced smaller declines or increases in LS and cluster coefficients.These webs also had larger proportions of intermediate and omnivore species.Webs with high PC-2 scores lost a greater proportion of species, showed an increase in C, and decreased path lengths.Webs with low PC-2 scores lost fewer species, experienced a decrease in C, and increased in path length.

DISCUSSION
Our simulations have illustrated that foundation species can play an important role in the assembly and collapse of food webs.By definition, foundation species influence community composition and functioning largely through non-trophic interactions (Ellison et al. 2005).Here, we have shown that the trophic position of the species that receive benefits (in this case a decrease in metabolic rate) from the presence of a foundation species can influence the food web assembly process and the response of a food web to the loss of a foundation species.When a foundation species lowered the metabolic rate of only basal species the resultant webs were complex and species-rich.In general, basal model webs also were robust to foundation species removals, retaining high species richness and complexity.On the other hand, when a foundation species lowered the metabolic rate of only consumer species (our consumer model), all species (total model), or no species (control model) the resultant webs were species-poor and the consumer webs had low complexity (i.e.low C, LS, clustering coefficient).Furthermore, the subsequent removal of the foundation species from the consumer, total, and control model webs resulted in a greater loss of species and complexity than in the basal model webs.
One potential explanation for the species-rich complex food webs produced by basal models and the species-poor simplified webs produced by the consumer and total models may be found in the population dynamics of the system.When a foundation species lowers the metabolic rate of the consumers (top predators and intermediate consumers in both the consumer and total models), consumer populations reach higher abundances, which in turn can lead to stronger predator-prey interactions (Holling 1965, Abrams andGinzberg 2000).Strong interactions can lead to unstable predator-prey dynamics and result in the extinction of both the predator and the prey species (May 1972, McCann et al. 1998).In the basal model, lower metabolic rates increased energy for growth and reproduction, allowing basal species to withstand transient dynamics of early assembly or low initial population abundances.Once gaining a foothold, even non-foundational basal species can provide multiple energy pathways to species at higher trophic levels.And once the foundation species was removed, the other basal species were already established and maintained energy pathways to higher trophic levels, limiting further extinctions.
This mechanism is also consistent with the standard facilitation model of succession (Connell and Slatyer 1977), where later-successional (facilitated) species can maintain high abundances even after early-successional species have disappeared.Two important differences, however, are that in the field, foundation species persist in the system much longer than early-successional species, and associated species composition changes dramatically following foundation species removal (e.g., Orwig et al. 2013).
In addition to the trophic position of the target species that a foundation species influences, the magnitude of the metabolic rates of the associated species in the absence of the foundation species (or more generally, the cost of not having the foundation species) was also important in determining food-web structure and the response of food webs to foundation species removal.When metabolic rates were highest in the absences of foundation species (the 10× treatment), webs lost the most species both during assembly and after removal of the foundation species.The 10× treatment also was the only one for which webs collapsed entirely (to zero species).This collapse was observed most frequently in the control webs, in which the foundation species did not have any non-trophic interactions with other species.Interestingly, basal model webs in the 10× metabolic rate group maintained species richness at levels similar to those seen in the lower metabolic rate treatments.This result is consistent with that seen in the food-web assembly dynamics, and implies that facilitation of basal species by foundation species can overcome even the highest metabolic rates (costs).Overall, our results suggest that foundation species that influence other basal species will result in robust food webs, whereas those that influence consumers lead to the loss of species and complexity both during the assembly process and after foundation species removal.Additionally, these effects are magnified when metabolic costs to other species increase in the absence of the foundation species.
In our models, foundation species exerted influence by lowering metabolic rates for certain species .This is only one type of non-trophic interaction that can occur in a food web, and it is likely that foundations species have many other non-trophic interactions and effects (e.g., providing refuge from predators, facilitating establishment; Kéfi et al. 2012) that deserve further exploration.In addition, in all of our models, foundation species had a positive influence on all species at similar trophic positions.In real food webs, however, this generalization is unlikely to hold, as foundation species can have different effects on species that share the same trophic position and may also have negative effects on some species in the food web (e.g., Ellison et al. 2005b, Sackett et al. 2010, Prevey et al. 2010, Kane et al. 2011).Furthermore, the effects of foundation species in our simulations are strongest when associated species do really poorly without the foundation species present (i.e., the 5× and10× metabolic treatments).This result implies that the role of a foundation species largely depends on the magnitude of its influence, but weak trophic (McKann et al. 1998, Neutel et al. 2002, Rooney and McCann 2012) and facilitative links (Allesina and Tang 2012) are also import in maintaining network structure and dynamics.Thus, measuring the influence of foundation species on other species in the food web through experimental removal studies (e.g., Ellison et al. 2010, Sackett et al. 2010) will continue to be an important component of understanding foundation species roles in the assembly and collapse of food webs.
Future exploration of foundation species in both modeled and real food webs should consider how foundations species differentially influence species in similar trophic positions, the threshold of metabolic rates (or other factors that foundations species influence) at which food webs respond, and non-trophic interactions that influence model parameters other than metabolic rate.Nonetheless, this first theoretical exploration of foundation species in a food-web context shows that we should look for foundation species to strongly influence basal species, leading to robust species-rich food webs that are the least susceptible to cascading extinctions when foundation species are lost.(blue), 5× (red), or 2× (orange) the baseline, allometrically-scaled metabolic rate (dashed line).
As the biomass of the foundation species increases, metabolic rate declines asymptotically to the baseline.These functions are the six metabolic rate treatments that we applied to the predatorprey model.

Fig. 1 .
Fig. 1.Schematic diagrams of the four foundation species models; A) control, B) basal, C) consumer, D) total.White nodes are basal foundation species, gray nodes are other basal species, and black nodes are consumers.Solid black lines with arrows represent trophic interactions and dashed lines are non-trophic interactions (i.e., reduction in metabolic rate).

Fig. 2 .
Fig.2.Saturating functions (eqn.6) relating metabolic rate to foundation species biomass.A) In the absence of a foundation species, species have the baseline, allometricallyscaled metabolic rate (dashed line; eqn.5b).Increasing the biomass of the foundation species results in an asymptotic decline in metabolic rate to 0.5× (green), 0.2× (magenta), or 0.01× (cyan) the baseline.B) When foundation species biomass = 0, species have metabolic rates 10×

Fig. 3 .
Fig. 3. ANCOVA plots illustrating species richness (A) and principal axis scores (B, C) of food-web structure after food-web assembly (at t = 2,000 modeled time steps) as a function of metabolic rate and the four types of foundation species models.Green lines and points correspond to the basal model, Pink = consumer model, Blue = total model, and Orange = control model.

Fig. 4 .
Fig. 4. Principal component biplots of food-web metrics for assembled food webs (at t = 2,000 modeled time steps).Illustrations along each PC axis depict representative individual webs.

Fig. 6 .
Fig. 6.Principal component biplots of standardized change in food-web metrics for food webs after foundation species removal (i.e., ∆z = z t=2001 -z t=4000 / z t=2001 ).Text along each PC axis show general change in food web complexity and richness associated with each axis.

Table 1 .
Model Variables

Table 2 .
Metrics of food-web structure

Table 3 .
Principal component loadings for food-web structure after food-web assembly (t = 2,000 modeled time steps) and after foundation species removal (t = 4,000 time steps).