Diverse interactions and ecosystem engineering stabilize community assembly

TL;DR

This study models ecological community assembly considering diverse species interactions and ecosystem engineering, revealing how engineers influence stability, diversity, and extinction dynamics in ecosystems.

Contribution

It introduces a novel assembly model incorporating ecosystem engineers, highlighting their nonlinear effects on community robustness and diversity.

Findings

Ecosystem engineers create niches that promote diversity.

Engineering effects can both stabilize and destabilize communities.

Redundancies in engineering effects lower colonization barriers.

Abstract

The complexity of an ecological community can be distilled into a network, where diverse interactions connect species in a web of dependencies. Species interact not only with each other but indirectly through environmental effects, however the role of these ecosystem engineers has not yet been considered in models of ecological networks. Here we explore the dynamics of ecosystem assembly, where the colonization and extinction of species within a community depends on the constraints imposed by trophic, service, and engineering dependencies. We show that our assembly model reproduces many key features of ecological systems, such as the role of generalists during assembly, realistic maximum trophic levels, and increased nestedness with higher frequencies of mutualisms. We find that ecosystem engineering has large and nonlinear effects on extinction rates, facilitating robustness by creating niche space, but at the same time increasing the magnitude of extinction cascades. We emphasize the importance of redundancies in engineered effects and show that such redundancy lowers the barriers to colonization, promoting community diversity. Together, our results suggest that ecological engineers may enhance community diversity while increasing persistence by facilitating colonization and limiting competitive exclusion.

Equations

$${\sigma }_i=c_{\mathrm{n}}{n}_i-c_{\mathrm{e}}{e}_i-c_{\mathrm{v}}{v}_i,$$

$$\mathrm{E}\{{\mathcal{M}}_{\mathbf{P}}\}={\mathcal{S}}_{\mathbf{P}}\eta \left(1-\frac{1}{e}\right),$$

$$\mathrm{E}\{p_{\mathrm{m}}\}=\frac{\eta }{{\mathcal{S}}_{\mathbf{P}}{\left(1+\eta -\frac{\eta }{e}\right)}^2}.$$

$$\mathrm{E}\{{\mathcal{M}}_{\mathbf{P}}{\}}_{\mathrm{redundant}}=\eta {\mathcal{S}}_{\mathbf{P}}\frac{e-2}{e},$$

$$\varphi =\frac{e-2}{e-1}\approx \mathrm{0.418.}$$

Full text

Significance Community assembly is constrained by interactions between and among species, many of which can have lasting effects on the environment. We explore the influence of these ecosystem engineers on colonization and extinction dynamics using a network model that includes trophic, mutualistic, and engineering dependencies between species and the abiotic environment. We find that ecosystem engineering can stabilize assembly particularly when multiple engineers have similar effects on the community.

To unravel nature’s secrets we must simplify its abundant complexities and idiosyncrasies. The layers of natural history giving rise to an ecological community can be distilled – among many forms – into a network, where nodes represent species and links represent interactions between them. Networks are generally constructed for one type of interaction, such as food webs capturing predation Paine1966 ; Dunne2002 ; Pascual2006 or pollination networks capturing a specific mutualistic interaction Bascompte2013 , and continues to lead to significant breakthroughs in our understanding of the dynamical consequences of community structure May1972 ; Gross2009 ; Allesina2012 , assembly Ponisio2017 , and coevolution Guimaraes2017 . Recent interest in ‘multilayer networks’ comprising multiple interaction types (multitype interactions) may provide additional insight into these processes Kefi2016 ; Pilosof2017 . However, interactions where species affect others by altering the environment in a lasting way have not yet been incorporated into models of ecological networks. These interactions, known as ecosystem engineering Lawton1994 ; OdlingSmee2013 or more generally niche construction OdlingSmee2013b ; Fukami2015 , are quite common in nature and exist in almost every ecosystem.

Diverse interactions occur not only between species but indirectly through the effects that species have on their environment Jones1994 ; Olff2009 ; OdlingSmee2013 . Elephants root out large saplings and small trees, enabling the formation and maintenance of grasslands Haynes2012 and creating habitat for smaller vertebrates Pringle2008 . Burrowing rodents create shelter and promote primary production by aerating the soil Reichman2002 , salmon and aquatic invertebrates create freshwater habitats by changing stream morphology Moore2006 , and leaf-cutter ants alter microclimates, influencing seedling survival and plant growth Meyer2011 . These examples illustrate ecosystem engineering, where the engineering organism alters the environment on timescales longer than its own Hastings2007 .

Ecosystem engineering not only impacts communities on ecological timescales, but has profoundly shaped the evolution of life on Earth Erwin2008 . For example, the emergence of multicellular cyanobacteria fundamentally altered the atmosphere during the Great Oxidation Event of the Proterozoic roughly 2.5 Byrs BP Erwin2008 ; Schirrmeister2013 , paving the way for the biological invasion of terrestrial habitats. In the oceans it is thought that rRNA and protein biogenesis of aquatic photoautotrophs drove the nitrogen:phosphorous ratio (the Redfield Ratio) to ca. 16:1 matching that of plankton Loladze2011 , illustrating that engineering clades can have much larger, sometimes global-scale effects.

The effect of the environment on species is commonly included in models of ecological dynamics Woodward2010 ; Brose2012 ; Gibert2019b due to its acknowledged importance and because it can – to first approximation – be easily systematized. By comparison the way in which species engineer the environment defies easy systemization due to the multitude of mechanisms by which engineering occurs. While interactions between species and the abiotic environment have been conceptually described Olff2009 ; Getz2011 , the absence of engineered effects in network models was addressed by Odling-Smee et al. OdlingSmee2013 , where they outlined a conceptual framework that included both species and abiotic compartments as nodes of a network, with links denoting both biotic and abiotic interactions.

Here we model the assembly of an ecological network where nodes represent ecological entities, including engineering species, non-engineering species, and the effects of the former on the environment, which we call abiotic modifiers. The links of the network that connect both species and modifiers represent trophic (eat interactions), service (need interactions), and engineering dependencies, respectively (Fig. 1; see Materials and Methods for a full description). Trophic interactions represent both predation as well as parasitism, whereas service interactions account for all non-trophic interactions such as pollination or seed dispersal. In our framework a traditional mutualism (such as a plant-pollinator interaction) consists of a service (need) interaction in one direction and a trophic (eat) interaction in the other. These multitype interactions between species and modifiers thus embed multiple dependent ecological sub-systems into a single network (Fig. 1). Modifiers in our framework overlap conceptually with the ‘abiotic compartments’ described in Odling-Smee et al. OdlingSmee2013 . Following Pillai et al. Pillai2011 , we do not track the abundances of biotic or abiotic entities but only track their presence or absence. We use this framework to explore the dynamics of ecosystem assembly, where the colonization and extinction of species within a community depends on the constraints imposed by the trophic, service, and engineering dependencies. We then show how observed network structures emerge from the process of assembly, compare their attributes with those from empirical systems, and examine the effects of ecosystem engineers.

Our results offer four key insights into the roles of multitype interactions and ecosystem engineering in driving community assembly. First, we show that the assembly of communities in the absence of engineering reproduces many features observed in empirical systems. These include changes in the proportion of generalists over the course of assembly that accord with measured data and trophic diversity similar to empirical observations. Second, we show that increasing the frequency of mutualistic interactions leads to the assembly of ecological networks that are more nested, a common feature of diverse mutualistic systems Bascompte2003 , but are also less robust. Our third key result shows that increasing the proportion of ecosystem engineers within a community has nonlinear effects on observed extinction rates. While we find that a low amount of engineering increases extinction rates, a high amount of engineering has the opposite effect. Finally we show that redundancies in engineered effects promote community diversity by lowering the barriers to colonization.

**Assembly without ecosystem engineering

**Communities assemble by random colonization from a source pool. A species from the source pool can colonize if it finds at least one resource that it can consume (one eat interaction is satisfied; cf. Ref. 33) and all of its non-trophic needs are met (all need interactions are satisfied). As such, the service interactions are assumed to be obligate, whereas trophic interactions are flexible. Following the establishment of an autotrophic base, the arrival of mixotrophs and lower trophic heterotrophs create opportunities for organisms occupying higher trophic levels to invade. This expanding niche space initially serves as an accelerator for community growth.

Following the initial colonization phase, extinctions begin to slow the rate of community growth. Primary extinctions occur by the competitive exclusion of species sharing similar resources. A species’ competition strength is determined by its interactions: competition strength is enhanced by the number of need interactions and penalized by its trophic generality (number of prey) and vulnerability (number of predators). Secondary extinctions occur when species lose its last trophic or any of its service requirements. See Fig. 1D, E for an illustration of the assembly process. As the colonization and extinction rates converge, the community reaches a steady state around which it oscillates (Fig. 2A). See Materials and Methods and Supporting Information, section I for a complete description of the assembly process.

Assembly of ecological communities in the absence of engineering results in interaction networks with structures consistent with empirical observations. As the community reaches steady state, we find that the connectance of trophic interactions ($C=L/S^{2}$, where $S$ is species richness and $L$ is the number of links) decays to a value similar to that of the source pool (Fig. S1). Decaying connectance has been documented in the assembly of mangrove communities Piechnik2008 , however this decay is a statistical inevitability, as early in the assembly process small food webs will have high link density, from which it can only decline. Compared to trophic networks constructed using the Niche model Williams2000 given similar species richness and connectance, our framework results in networks with degree distributions of similar means but with reduced variance (Supporting Information, section II).

Recent empirical work has suggested that generalist species may predominate early in assembly, whereas specialists colonize after a diverse resource base has accumulated Piechnik2008 ; Gravel2011 . Here the trophic generality of species $i$ is defined $G_{i}=k^{\rm in}_{i}/(L/S)$ where $k^{\rm in}_{i}$ is the in-degree, or number of resources consumed, by species $i$ Williams2000 . A species is classified as a generalist if $G_{i}>1$ and a specialist if $G_{i}<1$. If generalism is scaled to the steady state link density (see Supporting Information, section III), we observe that generalists dominate early in assembly, with an increase in specialists as assembly progresses (Fig. 2B). This confirms expectations from the trophic theory of island biogeography Gravel2011 , where early communities with lower richness are less likely to support specialist consumers than late, species-rich communities. At steady state the proportion of specialists is ca. 56%, similar to empirical observations of assembling food webs Piechnik2008 .

The role of specialists early in assembly is primarily due to the accumulation of autotrophic specialists. This is evident when we observe that the trophic level (TL) distribution early in assembly ($t=5$) has an average ${\rm TL}=1.6$. Four trophic levels are typically established by $t=50$, where colonization is still dominant, and by the time communities reach steady state the interaction networks are characterized by an average ${\rm TL_{max}}$ ($\pm$ standard deviation) $=11\pm 2.8$ (Fig. 2C). While the maximum trophic level is higher than that measured in most predator-prey systems Williams2002 , it is not unreasonable if parasitic interactions (which we do not differentiate from other consumers) are included Lafferty2006 . Overall, the most common trophic level among species at steady state is ca. ${\rm TL}=4.75$.

The distribution of trophic levels changes shape over the course of assembly. Early in assembly, we observe a skewed pyramidal structure, where most species feed from the base of the food web. At steady state, we observe that intermediate trophic levels dominate, with frequencies taking on an hourglass structure (purple bars, Fig. 2C). Compellingly, the trophic richness pyramids that we observe at steady state follow closely the hourglass distribution observed for empirical food webs and are less top-heavy than those produced by static food web models Turney2016 .

**Structure and dynamics of mutualisms

**Nested interactions, where specialist interactions are subsets of generalist interactions, are a distinguishing feature of mutualistic networks Bascompte2003 . A nested structure has been shown to maximize the structural stability of mutualistic networks Rohr2014 , emerge naturally via adaptive foraging behaviors Valdovinos2016 ; Valdovinos2019 and neutral processes Krishna2008 , and promote the influence of indirect effects in driving coevolutionary dynamics Guimaraes2017 . While models and experiments of trophic networks suggest that compartmentalization confers greater stabilizing properties Stouffer2011 ; Gilarranz2017 , interaction asymmetry among species may promote nestedness in both trophic Araujo2010 and mutualistic systems Pires2011 . Processes that operate on different temporal and spatial scales may have a significant influence on these observations Massol2011 . For example, over evolutionary time, coevolution and speciation may degrade nested structures in favor of modularity Ponisio2019 , and there is some evidence from Pleistocene food webs that geographic insularity may reinforce this process Yeakel2013 .

Does the assembly of ecological networks favor nestedness when mutualistic interactions are frequent? Increasing service dependencies (need interactions; see Fig. 1) promotes both service-resource and service-service dependencies. Consider how species with more service interactions compare to those with fewer. More service interactions i) increase a species’ competition strength, lowering its primary extinction risk while also ii) increasing inter-species dependencies and its secondary extinction risk. While mutualisms convey fitness advantages in order to evolve, the latter highlights the potential risk associated with losing mutualistic partners Bond1994 ; Colwell2012 . Indeed, the balance that mutualists must maintain with their partners may have large implications for the future of global biodiversity Dunn2009 .

We find that as we increase the frequency of mutualistic interactions, the assembled community at steady state becomes more nested (Fig. 3). In this case, nestedness emerges from the assembly process and provides structural robustness. The robustness can be observed by examining the exclusionary differences between species in a simple nested motif (Fig. 3, inset). In the trophic motif, species with high vulnerability (multiple predators) are at greater risk of primary extinction via competitive exclusion. This will result in the secondary extinction of the specialist consumer, rendering the nested structure prone to change. In our framework mutualistic networks are generally formed by composite interactions, where the consumer species is engaged in a trophic interaction while the resource species is engaged in a service interaction. As such, the consumer species becomes a trophic partner and the resource species gains the competitive advantages of the service. If the competitive advantages of services are greater than the costs of vulnerability (see Materials and Methods), it is the low vulnerability species with fewer trophic partners that is at greater risk of exclusion (Fig. 3, inset). Because its elimination will not cascade, the nested structure will be more resistant to change.

Our results also suggest that the addition of mutualistic interactions comes at a cost to the assembling community. Because mutualisms increase dependencies between species, they also increase the frequency of secondary extinctions (Fig. 3). Measuring persistence in terms of the proportion of time species are established in the network reveals that more frequent mutualisms leads on average to lower persistence. At the community-scale, lower average persistence implies greater species turnover. Observations of empirical systems appear to support our model predictions. For example, assembling plant-pollinator systems have demonstrated high rates of species and interaction turnover, both during the assembly process and at the steady state Ponisio2017 .

We emphasize that we have restricted ourselves to examining the effects of obligate mutualisms, although the importance of non-obligate mutualisms has long been recognized Ramos2012 ; Vieira2015 ; Valdovinos2016 ; Ponisio2017 ; Valdovinos2019 . While the inclusion of non-obligate mutualists will lower the likelihood of cascading effects in systems with higher frequencies of service interactions, the loss of obligate mutualistic partners will have larger dynamic consequences than the loss of more flexible non-obligate mutualistic partners. As such, we do not expect inclusion of non-obligate mutualisms to alter the qualitative nature of our findings.

**Assembly with ecosystem engineering

**The concept of ecosystem engineering, or more generally niche construction, has both encouraged an extended evolutionary synthesis Laland2015 while also garnering considerable controversy Gupta2017 ; Feldman2017 . Models that explore the effects of ecosystem engineering are relatively few, but have covered important ground Hastings2007 ; OdlingSmee2013 . For example, engineering has been shown to promote invasion Cuddington2004 , alter primary productivity Wright2004 , and change the selective environment over eco-evolutionary timescales Kylafis2008 ; Krakauer2009 which can lead to unexpected outcomes such as the fixation of deleterious alleles Laland1999 . On smaller scales, microbiota construct shared metabolitic resources that have a significant influence on microbial communities Kallus2017 , the dynamics of which may even serve as the missing ingredient stabilizing some complex ecological systems Muscarella2017 .

We next explore the effects of ecosystem engineering by allowing species to produce abiotic modifiers as additional nodes in the ecological network (Fig. 1). These modifier nodes produced by engineers can serve to fulfill resource or service requirements for other species. The parameter $\eta$ defines the mean number of modifiers produced per species, drawn from a Poisson distribution (see Materials and Methods for details). Increasing the frequency of engineering interactions both increases the number of engineering species (those species making $\geq 1$ modifier) and the number of modifiers per species. There are two characteristics of engineering that have particular relevance for community assembly: modifiers can linger in the community even after the species that produce them have been excluded, and more than one engineer can produce the same modifier such that engineering redundancies increase with $\eta$ (Supporting Information, section I).

Increasing engineering has significant consequences for community robustness, but these effects also are sensitive to the frequency of service interactions within the community. We measure community robustness by i) rates of primary versus secondary extinctions, ii) species persistence, and iii) steady state community diversity. All measures were averaged over each species within the community across assembly time.

As the number of engineers increase, mean rates of primary extinction are first elevated and then decline (Fig. 4A). This nonlinear effect of engineering on rates of primary extinction results from two competing forces. Increased production of abiotic modifiers supplies consumers with additional resources, limiting secondary extinctions and promoting persistence (Fig. 4B, C). However, the stabilization of consumers ultimately results in increased vulnerability of prey. Engineering dependencies are considered rare when the ratio of modifier nodes per species is $0<\eta\leq 0.5$. The cumulative effect in these species-rich/modifier-poor systems is increased competitive exclusion of prey and higher rates of primary extinction (Fig. 4A). Notably the presence of even a small number of engineers serves to limit the magnitude of secondary extinction cascades. Higher rates of primary extinction coupled with lower rates of secondary extinction mean that extinctions are common, but of limited magnitude such that disturbances are compartmentalized. As engineering becomes common ($\eta>0.5$) the available niche space expands, lowering competitive overlap and suppressing both primary and secondary extinctions.

Increasing the frequency of service interactions promotes service interactions between species and engineered modifiers (Fig. 1). A topical example of the latter is the habitat provided to invertebrates by the recently discovered rock-boring teredinid shipworm (Lithoredo abatanica) Shipway2019 . Here, freshwater invertebrates are serviced by the habitat modifications engineered by the shipworm, linking species indirectly via an abiotic effect (in our framework via a modifier node). As the frequency of service interactions increases, the negative effects associated with rare engineers is diminished (Fig. 4A). Increasing service interactions both elevates the competitive strength of species receiving services (from species and/or modifiers), while creating more interdependencies between and among species. As trophic interactions are replaced by service interactions, previously vulnerable species gain a competitive foothold and persist (Fig. 3, inset), lowering rates of primary extinctions (Fig. 4A). The costs of these added services to the community are an increased rate of secondary extinctions (Fig. 4B) and higher species turnover (Fig. 4C). Low rates of primary extinction coupled with high rates of secondary extinction mean that extinctions are less common but lead to larger cascades.

While the importance of engineering timescales has been emphasized previously Hastings2007 , redundant engineering has been assumed to be unimportant Lawton1994 . We argue that redundancy may be an important component of highly engineered systems, and particularly relevant when there exists a positive feedback between the effects of engineers on their fitness Cuddington2004 . The vast majority of contemporary ecosystem engineering case studies focus on single taxa, such that redundant engineers appear rare Lawton1994 . However if we consider longer timescales, increasing diversity of engineering clades may promote redundancy, and in some cases this may feed back to accelerate diversification OdlingSmee2013b . Such positive feedback mechanisms likely facilitated the global changes induced by cyanobacteria in the Proterozoic Erwin2008 ; Schirrmeister2013 among other large-scale engineering events in the history of life Erwin2008 . Engineering redundancies are likely important on shorter timescales as well. For example, diverse sessile epifauna on shelled gravels in shallow marine environments are facilitated by the engineering of their ancestors, such that the engineered effects of the clade determine the future fitness of descendants Kidwell1986 . In the microbiome, redundant engineering may be very common due to the influence of horizontal gene transfer in structuring metabolite production Polz2013 . In these systems, redundancy in the production of shared metabolitic resources may play a key role in community structure and dynamics Kallus2017 ; Muscarella2017 .

When there are few engineers, each modifier in the community tends to be unique to a particular engineering species. Engineering redundancies increase linearly with $\eta$ (Supporting Information, section I; Fig. S5), such that the loss of an engineer will not necessarily lead to the loss of engineered modifiers. We examine the effects of this redundancy by comparing our results to those produced by the same model, but where each modifier is uniquely produced by a single species. Surprisingly, the lack of engineering redundancies does not alter the general relationship between engineering and measures of community robustness (Fig. S6). However we find that redundancies play a central role in maintaining species diversity. When engineering redundancies are allowed, steady state community richness $S^{*}$ does not vary considerably with increasing service interactions and engineering (Fig. S7A). In contrast, when redundant engineering is not allowed, steady state community richness $S^{*}_{u}$ declines sharply (Figs. 4D, S7B).

Communities lacking redundancy have lower species richness because sparse interdependencies preclude colonization (Fig. S7C,D). Colonization occurs only when trophic and service dependencies are fulfilled. A species requiring multiple engineered modifiers, each uniquely produced, means that each required entity must precede colonization. This magnifies the role of priority effects in constraining assembly order Fukami2015 , precluding many species from colonizing. In contrast, redundancy increases the niche space available to species while minimizing priority effects by allowing multiple engineers to fulfill dependencies. Our results thus suggest that redundant engineers may play important roles in assembling ecosystems by lowering the barriers to colonization thereby promoting community diversity.

Together, the results of our model point to the importance of considering multitype interactions both between species and as mediated through changes to the environment via engineering. We suggest that including the effects of engineers, either explicitly as we have done here, or otherwise, is vital for understanding the inter-dependencies that define ecological systems. As past ecosystems have fundamentally altered the landscape on which contemporary communities interact, future ecosystems will be defined by the influence of engineering today. Understanding the role of ecosystem engineers is thus tantamount to understanding our own.

**Materials and Methods

** We model an ecological system with a network where nodes represent ecological entities such as populations of species and or the presence of abiotic modifiers affecting species such as (examples). Following Pilai et al. Pillai2011 , we do not track the abundances of entities but only track their presence or absence. The links of the network represent interactions between pairs of entities (x,y). We distinguish three types of such interactions: x eats y, x needs y to be present, x makes modifier y.

The assembly process entails two steps: first a source pool of species is created, followed by colonization/extinction into/from a local community. The model is initialized by creating $S$ species and $M=\eta S$ modifiers, such that $N=S+M$ is the average total number of entities and $\eta$ is the average number of modifiers per species in the system. For each pair of species (x,y) there is a probability $p_{e}$ that x eats y and probability $p_{n}$ that x needs y. For each pair of species x and modifier m, there is a probability $q_{e}$ that x eats m and a probability $q_{n}$ that species x needs modifier m. Additionally, each species makes a number of modifiers that is drawn from a Poisson distribution with mean $\mu=\eta e/(e-1)$ where $e$ is Euler’s number. Once the number of modifiers per species is determined, each modifier is assigned to a species independently. This means that multiple species may make the same modifier, and that there may be some modifiers that are not made by any species, which are eliminated from the pool.

In addition to interactions with ecosystem entities, there can be interactions with a basal resource, which is always present. The first species always eats this resource, such that there is always a primary producer in the pool. Other species eat the basal resource with probability $p_{e}$. Species with zero assigned trophic interactions are assumed to be primary producers. See Supporting Information, section I for additional details on defining the source pool.

We then consider the assembly of a community which at any time will contain a subset of entities in the pool and always the basal resource. In time, the entities in the community are updated following a set of rules. A species from the pool can colonize the community if the following conditions are met:

1. all entities that a species needs are present in the community, and
2. at least one entity that a species eats is present in the community. If a colonization event is possible, it occurs stochastically in time with rate $r_{\rm c}$.

An established species is at risk of extinction if it is not the strongest competitor at least one of its resources that it eats. We compute the competitive strength of species $i$ as

[TABLE]

where $n_{i}$ is the number of entities that species $i$ needs, $e_{i}$ is the number of entities from the pool that species $i$ can eat, and $v_{i}$ is the number of species in the community that eat species $i$. This captures the ecological intuition that mutualisms provide a fitness benefit, specialists are stronger competitors than generalists, and many predators entail an energetic cost. The coefficients $c_{\rm n},\leavevmode\nobreak\ c_{\rm e},\leavevmode\nobreak\ c_{\rm v}$ describe the relative effects of these contributions to competitive strength. In the following, we use the values $c_{\rm n}>c_{\rm e}>c_{\rm v}$, such that the competitive benefit of adding an additional mutualism is greater than the detriment incurred by adding another prey or predator. A species at risk of extinction leaves the community stochastically in time at rate $r_{e}$.

A modifier is present in the community whenever at least one species that makes the modifier is present. If a species that makes a modifier colonizes a community, the modifier is created immediately, however modifiers may persist for some time after the last species that makes the modifier goes extinct. Any modifier that has lost all of its makers disappears stochastically in time at rate $r_{m}$.

The model described here can be simulated efficiently with an event-driven simulation utilizing a Gillespie algorithm. In these types of simulations, one computes the rates $r_{j}$ of all possible events $j$ in a given step. One then selects the time at which the next event happens by drawing a random number from an exponential distribution with mean $1/\sum_{j}{r_{j}}$. At this time, an event occurs that is randomly selected from the set of possible events such that the probability of event $a$ is $r_{a}/\sum_{j}{r_{j}}$. The effect of the event is then realized and the list of possible events is updated for the next step. This algorithm is known to offer a much better approximation to the true stochastic continuous time process than a simulation in discrete time steps, while providing a much higher numerical efficiency Gillespie1977 . Simulations described in the main text have default parameterizations of $S=200$, $p_{\rm e}=0.01$, $c_{\rm n}=\pi$, $c_{\rm e}=\sqrt{2}$, $c_{\rm v}=1$, and $4000$ iterations.

Supplementary Information

Section I: Building the source pool

Here and henceforth, we refer to the assembly model presented in the main text as the ENIgMa model (E:eat, N:need, Ig:ignore, Ma:make). To initiate the ENIgMa assembly model, we must first constract the source pool, where each ecological entity (species and modifiers) is defined by its potential interactions with each other. The source pool interaction matrix $\bm{P}$ is generated by first setting the number of species in the pool $\mathcal{S}_{\bm{P}}$ and determining the number of modifiers $\mathcal{M}_{\bm{P}}$ that are made by ecosystem engineers. The resulting matrix is $\mathcal{N}_{\bm{P}}\times\mathcal{N}_{\bm{P}}$ where $\mathcal{N}_{\bm{P}}=\mathcal{S}_{\bm{P}}+\mathcal{M}_{\bm{P}}$, and is subdivided into four quadrants, only two of which play a role here: species-species interactions and species-modifier interactions (see Fig. 1). In these two quadrants, the expected frequency of eat interactions ${\rm E}\{p_{\rm e}\}$ and the expected frequency of need interactions ${\rm E}\{p_{\rm n}\}$ are free parameters, as is the expected number of modifiers made per species ${\rm E}\{\mathcal{M}_{i}\}=\eta$. Here and throughout, we simplify this parameter space by assuming that the frequency of eat and need interactions for species-species (SS) interactions and species-modifier (SM) interactions are equivalent, such that ${\rm E_{SS}}\{p_{\rm e}\}={\rm E_{SM}}\{p_{\rm e}\}$ and ${\rm E_{SS}}\{p_{\rm n}\}={\rm E_{SM}}\{p_{\rm n}\}$. For each species, a set number of modifiers is drawn from ${\rm Poiss}(\eta)$, such that the expected proportion of species that are engineers (species that make modifiers) is $1-{e}^{-\eta}$. If a particular modifier is randomly and independently drawn for a given engineer from a complete list of all possible modifiers, such that multiple species – with some probability – can make the same modifier, the expected number of modifiers is

[TABLE]

where $e$ is Euler’s number. The frequency of engineering (make) interactions is then calculated as

[TABLE]

Finally the frequency of the non-interaction is calculated as ${\rm E_{SS}}\{p_{\rm\varnothing}\}=1-{\rm E_{SS}}\{p_{\rm e}\}+{\rm E_{SS}}\{p_{\rm n}\}$ and ${\rm E_{SM}}\{p_{\rm\varnothing}\}=1-{\rm E_{SM}}\{p_{\rm e}\}+{\rm E_{SM}}\{p_{\rm n}\}+{\rm E_{SM}}\{p_{\rm m}\}$ for species-species and species-modifier interactions, respectively. Pairwise interactions are assigned randomly from these probabilities between species-species and species-modifiers independently in both quadrants, such that the source pool matrix has no imbued structure apart from the number of species, the number of modifiers, and the frequency of each directional interaction type. Each source pool is provided a basal resource (the first row/column). A species with a trophic interaction to this resource is identified as an autotroph (or mixotroph depending on its other trophic interactions). If they do not have service dependencies with other species/modifiers, it is these species that are uniquely able to initiate assembly.

We can determine analytically the expected number of unique versus redundant modifiers in the source pool. As the total number of modifiers is given in Eq. S1, the number of unique modifiers is given by ${\rm E}\{\mathcal{M}_{\bm{P}}\}_{{\rm unique}}=\mathcal{S}_{\bm{P}}\eta e^{-1}$. The number of redundant modifiers is then given as

[TABLE]

such that the proportion of redundant modifiers $\phi$ is

[TABLE]

Accordingly, we find that the number of redundant modifiers increases linearly with $\eta$, while the proportion of modifiers that are redundant is fixed. Figure S5A,B shows both analytical expectations and numerically-derived measures for ${\rm E}\{\mathcal{M}_{\bm{P}}\}_{{\rm redundant}}$ and $\phi$, respectively.

Section II: Comparison to niche model

We compared certain structural features of ENIgMa at steady state to those of the Niche Model Williams2000 . Comparisons were restricted to networks constructed in the absence of engineering because engineers introduce indirect effects that are not considered in static food web models, and may make such comparisons irrelevant. While there are many similarities, there are also some important differences, some of which are highlighted in the main text. While we consider a comparison of our framework with other food web models such as the Niche Model relevant, we emphasize that the motivations underlying both are distinct. Our approach is intended to provide a deeper understanding into how multitype dependencies between species and the environment impact the dynamics of community assembly. While capturing general qualitative features of empirical systems demonstrates that the dynamics we consider are ecologically relevant, the goal of our approach is distinct from that of static food web models, which aim to maximize structural similarities between model and empirical systems Williams2000 ; Williams2011 .

We compared steady state ecological networks that emerge from ENIgMa (described in Materials and Methods, main text) with food webs constructed from the Niche Model Williams2000 with similar species richness and connectance. Because species richness and connectance of the Niche Model are often altered by eliminating disconnected species, we compared i) species richness, ii) connectance, iii) mean species degree, iv) standard deviation of out-degree distributions, and v) standard deviation of in-degree distributions averaged across 1000 replicates for each model.

We found that all measures resulted in fairly similar values between ENIgMa and the Niche Model food webs with a some important differences (Figs. S2,S3). While similar, ENIgMa produces consistently lower values of connectance, mean species degree, as well as standard deviations of the in- and out-degree distributions. This means that the food webs produced by ENIgMa are more sparsely connected with less variance between species. These results were expected, as the Niche Model assumes systematically increasing dietary ranges with higher niche values, whereas the trophic interactions assigned to species in the source pool of ENIgMa are drawn independently. An important difference between the Niche Model and ENIgMa is that we do not distinguish between predators and parasites. A different framework known as the Inverse Niche Model Warren2010 has been proposed to address parasitic interactions. The Inverse Niche Model assumes increasing specialization with feeding hierarchies, which would serve to lower the average generality of species (lower degree). In addition, the Inverse Niche model outputs lower standard deviations of in- and out-degree distributions. Together these trends suggest that the qualitative structural differences that we observe for the assembly and Niche model may reflect an important structural distinction between food webs that do and do not include parasitic species.

Section III: Measures of generality

The trophic breadth of potential colonizers is thought to play an important role in community assembly. The definition of a specialist or generalist to some degree depends on the size and connectance of the larger food web. Trophic generality for a species $i$ is defined $G_{i}=k^{\rm in}_{i}/(L/S)$, where $k^{\rm in}_{i}$ is the in-degree, or number of resources consumed by species $i$ Williams2000 . A species is classified as a generalist if the number of its trophic interactions is greater than the average number of links per species, or $G_{i}>1$, and a specialist if $G_{i}<1$, where a community can be described by the proportion of specialists found therein. For interaction networks that are assembling over time, generality can be scaled by a number of different measures of $L/S$, and this has a large effect on our interpretation of the role of generality in community assembly. For instance, $L/S$ may be quantified by either including all autotrophic species or only autotrophic functional groups. Furthermore, the scaling of generality may be made with respect to the current state of the community at each point in time, or with respect to the community at steady state. For instance, in their investigation of assembling mangrove food webs, Piechnik et al. Piechnik2008 scaled trophic breadth to a standard steady state value of $L^{*}/S^{*}=0.2$ averaged across 102 food webs.

To examine how our assessment of the role of generalism over the course of assembly changes based on the application of different scalings, we employ three different measures of $L/S$ to calculate $G_{i}$:

1. $G_{i}^{{\rm all}}$, where $L$ accounts for all links in the food web and $S$ accounts for all species relative to each time interval in the assembly process (circles; Fig. S4b);
2. $G_{i}^{{\rm hetero}}$, where we consider only the links and species richness of heterotrophs, excluding autotrophs (points; Fig. S4b);
3. $G_{i}^{*}$, where $L$ and $S$ are measured with respect to the communities at steady state, which is most similar to the measure used to evaluate assembling mangrove food webs (diamonds; Fig. S4b). Whether trophic breadth is scaled to the current state of $L/S$ or the steady state value of $L^{*}/S^{*}$ has a large influence on the estimated proportion of generalists in the community, particularly when the size of the system is small. We observe that for $G_{i}^{{\rm all}}$, the system is initially assembled by specialist species, though over the course of assembly the proportion of specialists relative to generalists declines to intermediate values (circles representing the average over replicates in Fig. S4). If only the trophic links between non-autotrophs are considered as in $G_{i}^{{\rm hetero}}$, specialists still dominate early in assembly, but there is a greater range, such that some systems can be described by a mixed proportion of specialists and generalists (individual points representing independent replicates in Fig. S4).

The different normalizations by which generality is measured will impact the interpretation of both empirical and model systems alike. In our framework, species colonizing early in the assembly process are generalists compared to how the term is defined at the steady state, but they are functionally specialists with respect to the assembling community. For example, a species that is trophically connected to 10 resource species in the source pool may colonize a community where it is consuming a small subset of its potential range. As the community grows, that species may realize more of its trophic niche if those resource species subsequently colonize the system. To what end we label this species a generalist or specialist relative to the assembling community is thus subject to multiple interpretations.