Competition for nutrients increases invasion resistance during assembly of microbial communities.

The assembly of microbial communities through sequential invasions of microbial species is challenging to study experimentally. Here, I used genome-scale metabolic models of multiple species to model community assembly. Each such model represents all known biochemical reactions that a species uses to build biomass from nutrients in the environment. Species interactions in such models emerge from first biochemical principles, either through competition for environmental nutrients, or through cross-feeding on metabolic by-products excreted by resident species. I used these models to study 250 community assembly sequences. In each such sequence, a community changes through successive species invasions. During the 250 assembly sequences, communities become more species-rich and invasion-resistant. Resistance against both constructive and destructive invasions - those that entail species extinction - is associated with high community productivity, high biomass, and low concentrations of unused carbon. Competition for nutrients outweighs the influence of cross-feeding on the growth rate of individual species. In a community assembly network of all communities that arise during the 250 assembly sequences, some communities occur more often than expected by chance. These include invasion resistant "attractor" communities with high biomass that arise late in community assembly and persist preferentially because of their invasion resistance. Genome-scale metabolic models can reveal generic properties of microbial communities that are independent of the resident species and the environment.

INTRODUCTION

In his 1958 popular science book “The ecology of invasions by animals and plants” British ecologist Charles S. Elton set the foundation for invasion biology as an independent field of science (Elton, 1958). In this short book, Elton proposed a number of hypotheses that became the subject of intensive research in subsequent decades. Perhaps the most prominent of them is the “biotic resistance” or “diversity‐stability” hypothesis. It posits that more species‐rich communities are less susceptible to invasion by non‐native species. Pertinent empirical evidence accumulated in the decades since Elton's book was published. Some such evidence, including from grasslands, subtropical wetlands, riparian and agricultural ecosystems, supports the hypothesis (Boughton et al., 2011; Brown & Peet, 2003; Dukes, 2002; Kennedy et al., 2002; Naeem et al., 2000; Peltzer & MacLeod, 2014; Peng et al., 2019; Tilman, 1997). Other evidence argues that species‐rich communities are not less but more susceptible to invasion (Altieri et al., 2010; Brown & Peet, 2003; Fridley et al., 2004; Levine, 2000; Lonsdale, 1999; Stohlgren et al., 2003; Stohlgren, Barnett, et al., 2006; Stohlgren, Jarnevich, et al., 2006). To resolve this so‐called “invasion paradox”, some studies suggest that diversity helps reduce invasibility at small spatial scales, but increases it at large spatial scales (Brown & Peet, 2003; Fridley et al., 2004, 2007; Hui & Richardson, 2017; Richardson & Pysek, 2006). Here, I studied how species invasions help assemble microbial communities, where species invasions have been less well‐studied than in macroscopic organisms.

Several reasons are responsible for the relative neglect of invasion dynamics in microbial communities. First, such dynamics are difficult to observe experimentally through standard microbial ecology techniques, such as 16S rDNA sequencing, or bulk measurements of metabolite concentrations (Aguirre de Carcer, 2020; Nemergut et al., 2013). In addition, microbes can disperse globally and show complex mosaic phylogenies shaped by horizontal gene transfer (Nemergut et al., 2013; Soucy et al., 2015), which can blur the distinction between native and invading species. Moreover, even though microbial community ecology is a very active field, it has developed independently from invasion biology (Latombe et al., 2021). It asks different questions, such as about principles of species coexistence, convergent assembly trajectories, historical contingency (priority effects), and multistable equilibria (Bittleston et al., 2020; Blasche et al., 2021; Datta et al., 2016; Estrela et al., 2021; Friedman et al., 2017; Furman et al., 2020; Goldford et al., 2018; Gralka et al., 2020; Hoek et al., 2016; Jami et al., 2013; Kehe et al., 2021; Kelsic et al., 2015; Ratzke et al., 2020; Sanchez‐Gorostiaga et al., 2019; Wolfe et al., 2014).

Because it is challenging to observe species invasions and community assembly in microbes, mathematical and computational models remain important to identify general assembly principles. Such models have a long tradition in studies of community assembly (Drake, 1990; Drossel et al., 2001; Hui et al., 2021; Law & Morton, 1996; McKane, 2004; Minoarivelo & Hui, 2016, 2018; Post & Pimm, 1983). Pertinent models fall into broad classes, such as (generalized) Lotka‐Volterra models (Baiser et al., 2010; Gonze et al., 2018; Hofbauer & Sigmund, 1988; Kuntal et al., 2019; Minoarivelo & Hui, 2018; Mittelbach & McGill, 2019; Servan et al., 2018), consumer resource models (Cuddington & Hastings, 2016; Goldford et al., 2018; Marsland et al., 2019, 2020; Mittelbach & McGill, 2019; Valdovinos et al., 2018), and adaptive dynamics models (Hui et al., 2021).

Existing models have several limitations when applied to microbial community assembly. First, with notable exceptions (Brunner & Chia, 2019; Goldford et al., 2018; Marsland et al., 2019), they are designed to understand predator–prey or plant‐pollinator communities, but not microbial communities, which are governed by trophic interactions mediated by the abiotic environment (Bittleston et al., 2020). For example, microbial species frequently feed on small molecules that other microbial species excrete as by‐products of their metabolism. Second, existing models usually assume a constant environment that is not changed by community members. In contrast, microbial species can secrete multiple metabolic by‐products that change the chemical environment and thus construct new ecological niches for other microbes (Basan et al., 2015; Bittleston et al., 2020; Pacheco et al., 2019; San Roman & Wagner, 2018). Third, existing models generally allow only for pairwise interactions among species, whereas higher order interactions can be important (Grilli et al., 2017; Mickalide & Kuehn, 2019; Sanchez‐Gorostiaga et al., 2019). Finally and perhaps most importantly, existing models make specific and often ad hoc assumptions about which species interact. They are not suited to allow species interactions to emerge from first biological principles.

These limitations can be alleviated by the modelling framework used here (Khandelwal et al., 2013; Klitgord & Segre, 2010; Levy & Borenstein, 2013; Libby et al., 2019; Pacheco et al., 2019; San Roman & Wagner, 2021; Stolyar et al., 2007). It is built on genome‐scale metabolic models of microbial species (Amalric & Dehaene, 2019). Each such model represents all known metabolic biochemical reactions that take place in a given species. Its establishment requires a combination of genome sequence data and extensive biochemical information. Manual curation, together with increasingly advanced automatic model construction techniques have made such models available for hundreds of species (Devoid et al., 2013; Gu et al., 2019; Magnúsdóttir et al., 2017). With the computational method of flux balance analysis (FBA, see Methods), a genome‐scale model can be used to predict the instantaneous biomass growth rate of a species in a given chemical environment (Orth et al., 2010), as well as the amounts of metabolic by‐products that a microbe excretes into the environment. The predictions of this method have been experimentally validated in multiple experiments (Edwards et al., 2001; Feist et al., 2007; Segre et al., 2002; Varma & Palsson, 1994).

Genome‐scale metabolic modelling and FBA allow species interactions to emerge from basic biochemical principles. First, they allow species to compete by consuming the same nutrients available in the environment. Second, they can predict the identity and amount of metabolic by‐products that a species excretes and that can help other species grow. In doing so, they also allow a community's species to change their environment and create new ecological niches. Finally, they also allow higher order trophic species interactions to emerge naturally, because all species in a community have access to the same chemical environment.

Here, I use these methods to study the assembly of two different kinds of communities in silico. The first harbour a set of 100 human gut microbial species, and are assembled in an anaerobic chemical environment similar to a Western diet. The second harbour a subset of 100 random viable metabolisms (“species”), each of which contains a random complement of chemical reactions sampled from a known “universe” of such reactions realized in the biosphere, under the sole constraint that each such metabolism must be at least viable on glucose as the sole carbon source. I assemble the latter communities in an aerobic environment that contains only glucose as a carbon source. The two classes of metabolic systems are dramatically different, both in terms of their origins (biological/synthetic), and their environments (complex/simple, anaerobic/aerobic). Thus, any properties they share are probably generic properties of microbial communities assembled through trophic interactions, and not just peculiarities of a specific set of microbial species. For each type of community, I model 250 community assembly sequences. Each sequence starts out with a single species and proceeds through the sequential invasion of randomly chosen species from the species pool. I show that over time, both types of communities become more species‐rich, show higher biomass, and become more invasion resistant, largely because they become occupied by superior competitors. Some communities recur more often than expected by chance alone. They represent either recurrent “attractors” of community assembly or transient sources of new communities.

MATERIALS AND METHODS

Genome‐scale metabolic models and flux balance analysis

Here, each microbial species is represented by a genome‐scale metabolic model of bacterial metabolism (Becker et al., 2007; Feist et al., 2009; Schellenberger et al., 2010). Such a model comprises a set of stoichiometric equations for all known r (enzyme‐catalysed) chemical reactions that take place inside a living cell, and that involve the transformation of a total of m metabolites. These reactions also include a special class of so‐called exchange reactions by which metabolites such as nutrients are imported into a cell, and by which metabolic waste products are exported from a cell. Mathematically, a genome‐scale model consists of an stoichiometry matrix whose entries contains the stoichiometric coefficient of the i‐th metabolite in the j‐th reaction.

Genome‐scale models are assembled from a combination of biochemical information available through the literature and genome sequencing data. For a free‐living organism, they typically comprise more than 1000 reactions and metabolites. Each of the r reactions can proceed at some rate – the metabolic flux through reaction i – that has a lower bound and an upper bound , which are determined by factors such as the maximal rate at which the respective enzyme can catalyse the reaction, as well as by the thermodynamics of the reaction, such as whether the reaction is reversible or irreversible. Together, these fluxes define an ‐dimensional vector of metabolic fluxes.

If a genome‐scale model of a specific bacterial species exists, one can use the constraint‐based modelling technique flux balance analysis (FBA) to predict the biomass growth of the species in a given environment (Orth et al., 2010). To this end, FBA incorporates an (artificial) chemical reaction called the biomass growth reaction into a genome‐scale metabolic model. The substrates of this reaction are essential biomass precursors, such as amino acids, nucleotide building blocks, and lipids. Each of these substrates has a stoichiometric coefficient that reflects its proportion in a species' (experimentally determined) biomass composition. The biomass growth reaction then combines these substrates into biomass at biomass growth flux . More specifically, FBA assumes that a metabolism is in a steady state at which the concentrations of its internal metabolites do not change. This assumption would be met, for example, in a chemical environment like a chemostat that provides a continuous nutrient supply and continually removes metabolic waste products. FBA identifies a vector of metabolic fluxes that (i) conforms to the law of mass conservation in this steady‐state, and that (ii) maximizes the biomass growth flux in a given environment, subject to the constraints that exist for each metabolic flux . This task can be formulated as the dynamic programming problem (Orth et al., 2010). where the constraint reflects the steady‐state assumption. A solution to this optimization problem not only includes the maximally possible biomass growth rate, but also the rates at which exchange reactions must proceed to achieve this growth rate. These include the rates at which nutrients are consumed, and the rates at which by‐products of metabolisms are excreted into the environment.

In FBA, the environment is represented by the maximal rate at which different nutrients can be imported into a cell, that is, by the bounds on exchange reactions that represent the import of specific nutrients. Because it does not directly represent the concentrations of molecules and how they change when cells import nutrients and excrete waste products, FBA cannot be directly used to model the ecological dynamics of a microbial community over time. For this purpose, I use a simple and widely used extension of FBA called dynamic FBA (dFBA, Methods S1) (Chiu et al., 2014; Mahadevan et al., 2002; San Roman & Wagner, 2018; Varma & Palsson, 1994).

Community dynamics

I used dFBA to study the establishment of stable microbial communities through successful invasions from a pool of s species , each represented by a different genome‐scale metabolic model. To this end, I first chose a species at random (with equal probability for all species) from this pool, and introduced 0.01 g of this species' biomass into the chemostat. This choice ensures a sufficiently large number of initial cells to avoid the computationally costly complication of stochastic initial growth. Also, a substantially smaller inoculum would result in substantially longer amounts of time and computational cost for a community to reach equilibrium.

Subsequently, I simulated the population dynamics of this species with dFBA and determined every 24 h whether its biomass had either reached an equilibrium or whether the species had become extinct. (Extinction can take place if a species grows more slowly than the dilution rate of the chemostat.) In either case, I introduced 0.01 g of a new randomly chosen species (different from the resident species if that species had persisted) into the environment. Subsequently, I monitored the biomass growth of the resident species, and again tested for equilibrium every 24 h. Once such an equilibrium had been reached (or all species had gone extinct), I introduced another species, and so on.

I repeated this cycle of attempted invasions and community equilibration for 5000 h or for a maximum of 5 days of computation time on a single CPU, whichever was reached first (this 5 day maximum was reached in only 2% of community assemblies, which can happen when some communities take a very long time to reach equilibrium).

This procedure led to an assembly sequence of successive equilibrium microbial communities that hosted one or more species. A failed invasion led to the extinction of the invader. A successful invasion led to the establishment of the invader and either its coexistence with the resident community, or the extinction of one or more resident species. I recorded and analysed all equilibrium communities after a successful invasion, which are by definition communities whose species composition differed from the previous community in the assembly sequence.

Gut microbiome species

My starting point for the assembly of gut microbial communities is the assembly of gut organisms through reconstruction and analysis (AGORA) resource of metabolic models, which comprises 818 strains of gut bacteria (Magnúsdóttir et al., 2017). These models were reconstructed through a complex pipeline of genome annotation that is integrated with biochemical information and refined through comparative genomic data (Magnúsdóttir et al., 2017). FBA growth predictions of more than 50 of the models have been validated with experimental data (Magnúsdóttir et al., 2017; Noronha et al., 2019; Tramontano et al., 2018). Because information about species‐specific biomass composition is not available for most species, AGORA models use a general biomass growth reaction that yields biomass growth rates consistent with those observed in the human gut (Magnúsdóttir et al., 2017).

I downloaded the AGORA models from the Virtual Metabolic Human data base at https://www.vmh.life/#microbes/search on 9 April 2021 as “AGORA 1.03 without mucin” (Noronha et al., 2019). They come from 668 different microbial species and 214 different genera. To reduce the complexity of this data set while preserving as much of its diversity as possible, I sampled 100 species that can grow in the chemical environments I study from the complete set of models, such that each species comes from a different genus. Table S1 lists the resulting set of species that I used in this analysis. The modest sample size of 100 species also increases the likelihood that at least some species are sampled repeatedly during the assembly of different communities, which can help build informative assembly networks. The metabolic models for these species have 1286 ± 300 (1 SD) reactions, 871 ± 254 internal reactions, 415 ± 75 exchange reactions, and 1107 ± 192 metabolites.

I studied community assembly of these species on a Western diet (Magnúsdóttir et al., 2017), which comprises 164 nutrients, including many sugars, amino acids, fats, minerals, vitamins, and sources of fibre. Throughout, I assume anaerobic growth, to reflect the anaerobic or microaerobic nature of the human gut environment (Espey, 2013). The species I use (Table S1) excrete on average 12.7 ± 3.2 metabolites on a Western diet. They grow at a rate of μ ≈ 0.32 ± 0.2 h−1 and thus have a mean doubling time of h which is similar to doubling times observed in the mouse intestinal tract (Gibbons & Kapsimalis, 1967; Magnúsdóttir et al., 2017).

Continuous metabolic culture experiments that emulate microbial growth in the gut use a dilution rate D that is the inverse of the transit time of food through the human gut (Feria‐Gervasio et al., 2014; Macfarlane et al., 1998; Tottey et al., 2017). Because the median of this transit time is 60 h (Cummings et al., 1992), I chose a dilution rate of D = 1/60 ≈ 0.017 h−1. At this dilution rate, individual species show a mean growth rate of μ ≈ 0.32 ± 0.2 h−1, and 97 of the 100 species grow at a sufficient rate on a Western diet to persist on their own without being flushed out of the model gut.

Random viable metabolisms

To complement my analysis of communities assembled from human gut microbes, I also studied random viable metabolisms. Each such metabolism comprises a set of biochemical reactions that can produce all of a cell's essential biomass precursors and thus support biomass growth in a given environment, but that contain an otherwise random set of biochemical reactions samples from a “pan‐metabolism” or “universe” of such reactions that exists in the biosphere (Methods S1). The environment I use here is a minimal aerobic environment in which glucose is the only available carbon source.

RESULTS

Species richness increases but unevenness does not during community assembly

Following previous experimental work to establish a gut microbiome in vitro, I model the gut as a continuous microbial culture, that is, a chemostat to which nutrients are continually added and from which waste is continually removed (Feria‐Gervasio et al., 2014; Macfarlane et al., 1998; Tottey et al., 2017). Based on the experimentally observed residence time of food in the gut, I use a dilution rate of D = 0.017 h−1 for this purpose. I assume that the nutrient composition is similar to a previously described Western diet, which comprises 164 different nutrients (Methods S1). I assemble communities from a “regional” pool of 100 microbial species, each of which is represented by a genome‐scale metabolic model that comprises on average 1286 ± 300 (1 Sd) reactions that interconvert 1107 ± 192 metabolites.

I begin the assembly process by choosing a species from the pool at random, and introduce a small amount of biomass (0.01 g) of the species into the chemostat. Subsequently, I simulate the population dynamics of this species with dynamic FBA (Methods S1). This extension of FBA permits simulation of not only the biomass growth rate of one or more resident species, but also their consumption rate of each nutrient, and their excretion rate of each metabolic byproduct. I determine every 24 h whether the species' biomass has reached an equilibrium or whether the species had become extinct. In either case, I introduce 0.01 g of a new randomly chosen species, monitor the biomass of growth of the (one or two) resident species, and again test for equilibrium every 24 h (Methods S1). Such an equilibrium may again involve extinction of all, one, or none of the resident species. If the invading species is part of the new equilibrium community, I call its invasion successful. Once such an equilibrium has been reached, I introduce a third species, and so on, repeating the cycle of attempted invasions and community equilibration for a maximum of 5000 h (Methods S1).

This process leads to an assembly sequence of successive equilibrium microbial communities that host one or more species. I record for further analysis only equilibrium communities whose species composition is different from that of the previous communities, that is, communities in which the most recent invasion has been successful. For brevity, I will refer to them simply as communities. Figure 1a shows an example of the first 3000 h of such an assembly sequence, involving eight invasions numbered by integers. Invasions 1 and 2 establish a two‐species community. However, when a third species invades at approximately 100 h, it drives the two resident species to extinction while becoming established. I refer to such a successful invasion, which leads to the extinction of one or more resident species as a destructive invasion. At approximately 1000 h, a fourth species invades, which then coexists with the resident species. I refer to such a successful invasion, which leads to no extinction, as a constructive invasion. The example also illustrates that later communities need not necessarily have more species than earlier communities, because some species may go extinct after a successful invasion, for example, as a result of the third invasion in Figure 1a.

FIGURE 1

Assembly sequences create communities with multiple species. (a) Example of the first 3000 h (inset: 1000 h) of an assembly sequence involving the eight species on the left. The vertical axis shows the biomass [g] of each of species. Numbers inside the plot indicate the eight invasions that occur during this time period. Invasions 1–4, 6, and 7 are successful, whereas invasions 5 and 8 fail. Invasions 1, 2, 4, and 6 are constructive, whereas invasions 3 and 7 are destructive. Species names are as reported in (Magnúsdóttir et al., 2017; Noronha et al., 2019). (b) Distribution of the number of equilibrium communities that form during the 250 assembly sequences I simulate. (c) Mean number of species per community, (d) richness‐adjusted Shannon diversity, and (e) Berger‐Parker dominance (vertical axes, Methods S1), when communities are ranked by their order of appearance (first, second, etc, horizontal axes) in an assembly sequence. Note that Shannon diversity in (d) can only be calculated for communities of two or more species. Error bars in panels (c–e) indicate one standard error of the mean, and are often too small to be visible

Overall, I created 250 such assembly sequences, each of which lasted for a maximum of 5000 h (Methods S1). Some assembly sequences lead to many more communities than others, with a mean number of 9.5 ± 2.8 communities per assembly sequence (Figure 1b). Each of these communities can be ranked according to the order in which it appears in the sequence. For most of my analyses, I will pool the first, second, etc. community from different assembly sequences to analyse their properties jointly. In other words, I will quantify time during an assembly sequence not in hours since the beginning of the sequence, but through the first, second, third, etc community in the assembly sequence. To obtain statistically meaningful results below, I will analyse only the first 13 communities produced by my assembly sequences, because fewer than 30 assembly sequences generated 13 or more communities (Figure 1b).

Although the species richness – the number S of species per community – need not necessarily increase in consecutive communities, it increases on average when studied over multiple assembly sequences (Figure 1c). That is, whereas the first community in an assembly sequence comprises only one species, the last (13th) community comprises 6.9 ± 1.6 species, and the mean number of species per community steadily increases between these values. The evenness of the species abundance does not decrease substantially during community assembly (Figure 1d), and the dominant species accounts for an ever‐decreasing fraction of biomass (Figure 1e).

Communities become more resilient to invasions over time

Resilience is a system's ability to persist in the face of perturbations (Holling, 1973). The perturbations I study are species invasions, and I ask how likely it is that a community's species composition remains unchanged when it is invaded by a new species. To this end, I first quantify the number of unsuccessful invasions between two successful invasions, that is, the number of invasions that are “repelled” by a community before another community emerges. Figure 2a shows this number, averaged over the 250 assembly sequences I study. The establishment of the second community is associated with 0.11 failed invasions, a number that increases until it reaches a maximum at 0.88 failed invasions to establish the sixth community, and levels off thereafter. By this measure, invasion resistance increases by a factor eight during community assembly. In addition, I also evaluate the time that elapses between successive communities, which also increases substantially (Figure 2b). Specifically, while on average 35 min elapse to get the first (single‐species) community established, 614 min elapse on average between the sixth and the seventh community, a greater than 17‐fold increase. After reaching a maximum at the seventh community, this transition time decreases modestly, but remains substantially above that for early communities. (Its modest decline is caused predominantly by destructive invasions; Results S1, Figure S3). Although later communities in an assembly sequence take a longer time to reach equilibrium (Spearman's r = 0.32, p < 3 × 10−58; n = 2367), this association cannot solely explain the overall time increase between successful invasions. The reason is that the time between successive communities still increases when controlling for this equilibration time (partial Spearman's r = 0.21, p = 9.3 × 10−26; n = 2367; controlled for equilibration time). Not surprisingly, the greater the number of unsuccessful invasions, the greater is also the time needed between two successive communities (Spearman's r = 0.57; p = 6.9 × 10−203; n = 2367), and this association also persists when controlling for time to equilibrium (partial Spearman's r = 0.93, p < 10−200; n = 2367).

FIGURE 2

Invasion resistance increases early during community assembly. Vertical axes show (a) the mean number of failed invasions, and (b) the mean time between successive equilibrium communities during assembly of gut microbial communities), when communities are ranked by their order of appearance (first, second, etc) in an assembly sequence (horizontal axes). (c and d) like (a) and (b), respectively, but for random viable metabolisms. Data are reported up to a number of communities generated by at least 30 of 250 assembly sequences. This number is greater for gut microbial communities (13) than for random viable metabolisms (11). Error bars indicate one standard error of the mean

Like all organisms, gut microbes have a unique evolutionary and ecological history, and this history may influence the properties of their communities. My observations thus far may be a byproduct of this history, or they may be generic properties of microbial communities defined by their metabolic abilities. To find out, I turned to a complementary modelling approach that explores the assembly of communities from “species” that are defined by random viable metabolisms. Each such metabolism comprises a random subset of chemical reactions from a much larger reaction “universe” or “pan‐metabolic” set of reactions that are realized in some organism of the biosphere, subject to the requirement that the subset must be able to sustain life in a specific well‐defined chemical environment. Specifically, it must be able to synthesize all essential biomass precursors of a bacterial cell from the nutrients provided in this environment.

I use previously described procedures to create 100 such random viable metabolisms (Methods S1), and consider each such metabolism as a “species” for the purpose of community assembly. I created each of these species such that it has a number of internal chemical reactions that is identical to the average number of reactions of the 100 gut microbial species (Methods S1). Although communities assembled from such random viable metabolisms are not realized in the biosphere, they can help find out whether any community property may be a generic property of organisms interacting metabolically. In addition, they allow one to circumvent the highly complex nutrient requirements of gut microbes (Tramontano et al., 2018), because random viable metabolisms can be designed to be viable in specific nutrient environments. I created the 100 random viable metabolisms to be viable in an aerobic glucose minimal environment where glucose is the only source of carbon and energy (Methods S1).

With a procedure identical to that I used for the gut microbes, I then created 250 community assembly sequences that lasted for up to 5000 h. Each assembly sequence led to a number of communities that varied broadly around an average of 8.4 ± 2.6 communities (Figure S1A). This number is somewhat lower than the 9.5 ± 2.8 communities per assembly sequence for gut microbes. The average number of species in a community increased, albeit to a lower number of 4.8 ± 1.9 species in the last community (Figure S1B). The fewer communities and smaller number of species are not surprising, because they assembled in a simpler environment that can sustain fewer species. Similar to the gut microbiota, the unevenness of the species distributions did not increase dramatically over time (Figure S1C, D).

More importantly, both measures of invasion resistance increase dramatically early during an assembly sequence. Specifically, the number of failed invasions increases from an average of 0.25 between the first and second community by a factor 4.5 to an average of 1.12 between the fourth and fifth community, and stays approximately constant thereafter (Figure 2c). The time elapsed between consecutive communities also increases rapidly and early during an assembly sequence from 24.7 min to 732.5, that is, by a factor 29.7. It then declines modestly for late communities, but still remains far above the time needed to establish the first community. Once again, the two quantifiers are highly positively associated (Spearman's r = 0.49, p < 9 × 10−102; n = 2090), and even more so when one controls for the time that individual communities need to reach equilibrium (partial Spearman's r = 0.88, p < 10−100; n = 2090). Taken together, these observations show that properties such as the increase in invasion resistance during early community assembly are generic properties of microbes that interact metabolically through environmental nutrients and metabolic by‐products. They are also insensitive to the environment.

Productivity is associated with resistance to destructive invasions

I next studied the likely causes of a community's resistance to invasions. I first focused on resistance to destructive invasions, which I quantified as the number of failed invasions a community experiences before a destructive invasion occurs. I hypothesized that a community's productivity is important for such resistance. If an invader grows more rapidly than some members of the community, either on nutrients that exist in the environment or on nutrients that other species excrete as waste products, then it may cause a destructive invasion in which these members go extinct. As time progresses during an assembly sequence, ever‐superior competitors may replace inferior resident species in this way. The incidence of such destructive invasions may also depend on the earliest stage of an assembly sequence. Because all species are drawn from a finite species pool with a given distribution of biomass growth rates (Table S1), if the first invading species happens to grow slowly, it is more likely to be outcompeted by the next invading species. Over time, successful destructive invasions may become more rare, because the resident species increasingly consist of strong competitors.

This hypothesis creates several predictions. First, for it to be true, the communities I assemble must become increasingly more productive. This is indeed the case, because the equilibrium biomass of successive communities increases during an assembly sequence, both for gut microbes (Figure 3a) and for random viable networks (Figure S4A). Second, the hypothesis entails that communities become more invasion resistant as their biomass increases. This is also the case, and for both measures of invasion resistance, that is, for the number of failed invasions (Spearman's r = 0.28, p = 4.8 × 10−12; n = 604) and the time between successive communities (partial Spearman's r = 0.19, p = 2 × 10−6; n = 604, controlled for time to equilibrium).

FIGURE 3

Community productivity over time and productivity changes caused by destructive and constructive invasions. (a) Mean community biomass, (b) biomass growth rate of the dominant species on a Western diet, (c) mean number of metabolites in the environment, and (d) mean concentration [mmol] of one carbon (C1) units are shown on the vertical axes, and plotted against communities ranked by their order of appearance (first, second, etc) in an assembly sequence on the horizontal axis. Data are reported up to a number of communities generated by at least 30 of 250 assembly sequences. Error bars indicate one standard error of the mean, and are often too small to be visible. (e) Box‐whisker plot of the biomass change caused by constructive and destructive invasions; (f) box‐whisker plot of the reduction in C1 concentration caused by constructive and destructive invasions. Horizontal lines inside each box correspond to the median, boxes extends from the first to the third quartile of the data, and whiskers extend from the box through 1.5 times the interquartile range

A third prediction is that this increase in productivity also characterizes the dominant members of each community. This is indeed the case, both for communities of gut microbes (Figure 3b) and for random viable metabolisms (Figure S4B). In addition, with increasing productivity of this dominant member, communities become more invasion resistant (number of failed invasions: Spearman's r = 0.22, p = 3.8 × 10−8; n = 604; and the time between successive communities (partial Spearman's r = 0.13, p = 10−3; n = 604, controlled for time to equilibrium).

A fourth prediction relates to the first species that becomes established during an assembly sequence. If my hypothesis is correct, then the time needed to establish a second community through a destructive invasion should be smaller if the first invader grows slowly than if it grows rapidly. In other words, we would expect a positive correlation between the growth rate of the first invader and this time. This correlation indeed exists (Spearman's r = 0.83, p < 1.2 × 10−4, n = 15). Taken together, these observations suggest that more productive communities are more resistant to destructive invasions. Thus, increased community productivity (Figures 3a, b) is important to help explain increasing invasion resistance.

Decreasing carbon concentration is associated with resistance to constructive invasions

I next turned to examine resistance to constructive invasions, that is, the number of failed invasions a community experiences before a successful constructive invasion occurs. Such invasions can have at least two causes. First, in an environment providing multiple nutrients, an invader may exploit a resource that is not yet fully utilized by the resident species, that is, it may invade an unoccupied ecological niche. Second, the invader may become established by feeding on metabolic waste products of other species. Both causes may operate simultaneously.

To study any changes in resistance to constructive invasions, it is relevant that individual gut microbial species, when growing by themselves on a Western diet, excrete multiple metabolites (12.7 ± 3.2 metabolites per species). Not surprisingly then, during community assembly, the nutrient environment becomes chemically more and more complex, as communities harbour more and more species (Figure 1c), and these species excrete more and more metabolites. Specifically, during an assembly sequence the average number of metabolites in the environment increases from 126.1 ± 11.2 to 144.3 ± 13.8 (Figure 3c). However, the total amount of carbon in the environment actually declines (Figure 3d). To compute this amount, I determined the environmental concentration of each carbon‐containing nutrient with a known stoichiometry, and multiplied it with the number of carbon atoms the molecule contains. This procedure takes into account that some molecules contain much more carbon than others.

I hypothesized that environments with less carbon provide fewer opportunities for constructive invasions and would thus be more invasion resistant. If so, then the carbon content of the environment should be negatively correlated with invasion resistance. This is indeed the case for both measures of invasion resistance, that is, for the number of failed invasions (Spearman's r = −0.35, p = 2.8 × 10−44; n = 1513), and for the time between successive communities (partial Spearman's r = −0.34, p = 2.2 × 10−43; n = 1513, controlled for time to equilibrium).

In sum, when the resistance of communities to constructive and destructive invasions increases over time during an assembly sequence, it does so for reasons related to productivity and efficient resource use. A more productive late community is less easy to invade by a superior competitor that causes resident species extinctions. By the same token, such a community uses nutrients more efficiently, leaves fewer unused nutrients in the environment, and thus reduces opportunities for constructive invasions. Both destructive and constructive invasions tend to increase productivity and resource efficiency, but they do so to a different extent. Specifically, destructive invasions increase equilibrium community biomass by twice as much as constructive invasions (by 0.34 vs. 0.17 g, Figure 3e), a difference that is highly significant (p = 9.3 × 10−5, two‐tailed Mann–Whitney U test). Conversely, destructive and constructive invasions reduce the amount of free carbon by an average of 5.6 and 3.7 mmol, respectively, a difference that is also highly significant (Figure 3f , p = 1.5 × 10−5, two‐tailed Mann–Whitney U test). Thus, although fewer destructive invasions occur during community assembly (Figure S2A), they affect community productivity more drastically. Similar observations hold for communities assembled from random viable metabolisms, where constructive invasions also outnumber destructive invasions (Results S2, Figures S2B and S4).

The importance of competition is underscored by a complementary analysis that asks whether a species' biomass growth is affected predominantly negatively (competition) or positively (facilitation) by other members of the same community. It shows that at least 80% of pairwise interactions (and even more higher‐order interactions) are competitive, both for gut microbes and randomly viable metabolisms (Figures S5 and S6, Results S3). At the same time, cross‐feeding occurs between at least 70% of species pairs, but its effect on biomass growth is not sufficient to outweigh competition (Figure S7, Results S3). Even so, cross‐feeding contributes to community diversity, especially for random viable metabolisms. The reason is that in a glucose‐minimal environment, cross‐feeding is essential to sustain multi‐species communities.

A community assembly network helps identify highly invasion resistant communities as attractors of community assembly

Even species drawn from a modest pool of 100 species can potentially form 2100 = 1.3 × 1030 communities, many more than can be explored. Even if I also consider all unsuccessful invasions (which lead to communities where not all resident species can coexist), the 250 assembly sequencies I studied sample merely a fraction 3408/1.3 × 1030 = 2.6 × 10−27 of all possible communities. Despite this highly sparse sampling of all possible communities, some communities of two or more species recur repeatedly across different assembly sequences (Figure 4). To study such recurrences systematically, I constructed a directed community assembly network (graph) by pooling information from all 250 assembly sequences (Figure 4a,b). Each node in this network is a community of coexisting species, and two communities C and C are connected by a directed edge if community C emerged from community C via a successful species invasion during one of the assembly sequences. If this invasion was constructive, then C contains all species of C and additionally the invading species. If the invasion was destructive, C contains at most as many species as C , and possibly fewer. The resulting network contains N = 2123 nodes (communities) and K = 2112 directed edges. Its most basic property is its degree distribution, that is, the number k of communities that are neighbours of any one community, that is, connected by an edge to it (Figures 4c–f). Because the assembly network is directed, one needs to distinguish between the indegree k in and the outdegree k out of each community (k = k in + k out), that is, the number of edges emanating and terminating at any one community.

FIGURE 4

A community assembly graph shows that some communities occur repeatedly across assembly sequences. (a) Largest (“giant”) component (N = 1117 nodes, K = 1139 edges) of the community assembly graph, which comprises N = 2123 nodes (communities) and K = 2112 edges, and fragments into 50 weakly connected components. Each circle corresponds to one equilibrium community, and two communities C and C are connected by a directed edge if community C emerged from community C via a successful species invasion during an assembly sequences. Circle size corresponds to community biomass, and darker shading indicates higher indegree. (b) Magnification of the boxed area from (a). Note that several high biomass nodes have darker shading, which indicates that community assembly reached such nodes more than once. The bidirectional arrows indicate communities that can be reached from each other via a successful species invasion. (c) Distribution of indegrees for all communities. (d) Top 10 communities of two or more species ranked by indegree, together with the probability (p‐value) of observing their indegree by chance alone and the species identifiers (Table S1) of their constituent species. (e) Distribution of outdegree for all communities. (f) Top 10 communities of two or more species ranked by outdegree, together with the probability (p‐value) of observing their outdegree by chance alone and the species identifiers (Table S1) of their constituent species. Note that all listed communities consist only of a single species

For communities (nodes) of two or more species, the indegree quantifies the number of times the community occurs in the 250 assembly sequences. To ask whether the recurrence of some communities could be expected by chance alone, I employed a simple statistical test which assumes that all species can form persistent communities (Methods S1). For example, a five‐species community of species 1, 2 18, 28, and 32 (Figure 4d, Table S1 for species identifiers) occurred twice. Only 341 of the 2123 communities I sampled comprise five species, which is a vanishing fraction of the 7.5 × 107 possible five species communities. Not surprisingly then, the likelihood of encountering a community of five species twice is not expected by chance alone (p = 5.5 × 10−12, Figure 4d, Methods S1).

Several other communities also recur more frequently than expected by chance alone (Figure 4d). This suggests that some biological properties of these communities are responsible for this recurrence. To find such properties, I studied the subset of communities that vary in their degree, that is, where k in > 1 or k out > 1. I found that, the number of times that a community recurs (k in) is positively associated with its biomass (Spearman's r = 0.79, p = 2.4 × 10−29, n = 131), the biomass of its dominant species (Spearman's r = 0.63, p = 3.9 × 10−16, n = 131), and the biomass growth rate of this species on a Western diet (Spearman's r = 0.67, p = 4.1 × 10−18, n = 131). Thus, recurrent communities are productive communities. They also tend to leave lower amounts of carbon (C1 units) in the environment that could be used for biomass growth (Spearman's r = −0.76, p = 2.5 × 10−26, n = 131).

Perhaps surprisingly, a community's recurrence is positively associated with the average number of failed invasions that were needed to establish it from its predecessor (Spearman's r = 0.52, p = 1.96 × 10−10, n = 131), and with the average time needed to establish them from the communities preceding them (Spearman's r = 0.79, p = 5.3 × 10−29, n = 131). From this point of view, recurrent communities are more difficult to reach from other communities. This is, however, consistent with the fact that they are highly productive, because highly productive communities arise late during an assembly sequence, when the invasion resistance of communities has also become high. Indeed, recurrent communities do tend to arise late during an assembly sequence (Spearman's r = 0.92, p = 2.9 × 10−54, n = 131). In addition, once established, recurrent communities are less likely to be replaced by other communities, that is, their invasion resistance is higher (Spearman's r = 0.53 and 0.48 between k in and number of failed invasions until the next community, and time to next community, respectively, p < 5.1 × 10−9, n = 131).

As opposed to the indegree, the outdegree of a community in the assembly network indicates how often the community gives rise to other communities through successful invasions (Figures 4e, f, S8). It can thus also be viewed as a measure of a community's transience. The communities with the highest outdegrees all consist of single species. For example, species 63 has an outdegree of k out = 9 (Figure 4f), that is, it gives rise to nine other communities. This high outdegree is not expected by chance alone (p = 1.1 × 10−6, Methods S1). Multiple other communities have high outdegrees, which raises again a question about biological correlates of community transience. Communities with high outdegree tend to have low biomass (Spearman's r = −0.35, p = 4.3 × 10−5, n = 131), and low biomass growth rate of the dominant species on a Western diet (Spearman's r = −0.19, p = 0.03, n = 131). They also tend to leave high amounts of carbon (C1) units unused that could be used for the biomass growth of invading species (Spearman's r = 0.38, p = 7 × 10−6, n = 131). In other words, these communities have low productivity. Not surprisingly then, they tend to occur early during an assembly sequence (Spearman's r = −0.42 between k out and a community's position in an assembly sequence; p = 5.7 × 10−7, n = 131). In addition, communities with high outdegree are easily reached from their predecessors. The average number of failed invasions leading to them tends to be small (Spearman's r = −0.41, p = 8.7 × 10−7, n = 131), as does the average time needed for transitions between their predecessor and themselves (Spearman's r = −0.22, p = 0.01, n = 131). Conversely, few failed invasions and little time are on average needed to create their successor communities (Spearman's r = −0.21 between k out and number of failed invasions to successor, p = 0.015, r = −0.23 between k out and time to successor, p = 0.089; n = 131). In these analyses, neither k in nor k out were significantly associated with several other properties that I examined in previous sections, that is, the incidence of competitive, neutral, facilitative, and cross‐feeding interactions. Analogous associations exist for random viable metabolisms (Results S4).

In summary, communities that occur recurrently during community assembly are highly productive, occur late during an assembly sequence, are more difficult to reach, but are also more invasion resistant. Conversely, transient communities are less productive, occur earlier during an assembly sequence, are easier to reach, but are also less invasion resistant. Transient communities are in several ways the opposite of recurrent communities. Indeed, the indegree and the outdegree of communities are negatively associated (Spearman's r = −0.31 p = 2.9 × 10−4, n = 131).

DISCUSSION

For both gut microbes and random viable “species”, I find that communities become more invasion resistant, especially during community assembly. These observations are consistent with earlier modelling work (Drake, 1990; Law & Morton, 1996; Post & Pimm, 1983). However, in contrast to the Lotka‐Volterra models of earlier work, species interactions in my genome‐scale metabolic modelling framework emerge from first biochemical principles. My observations thus suggest that increased invasion resistance is a generic feature of microbial communities in which trophic interactions mediated by environmental nutrients are prevalent.

The models I use also help explain why invasion resistance increases. At the center of this explanation is a community's biomass productivity. It helps explain resistance against both constructive and destructive invasions, but for different reasons. Productive communities harbour few opportunities for constructive invasions, because they harbour few resources that are underused by resident species, which an invader can exploit without interfering with the persistence of a resident. Conversely, productive communities also harbour few opportunities for destructive invasions, because one or more of their resident species convert nutrients efficiently into biomass, such that few invaders can outcompete them.

In the communities I study, most pairwise species interactions are competitive (>80%) rather than facilitative (<17%), that is, a species' growth rate is reduced rather than increased by the presence of another species. This is not surprising, because all species can survive on the externally provided nutrients, and are thus competing for these nutrients. It holds even though 70% or more of species pairs cross‐feed. The prevalence of competition together with extensive cross‐feeding imply that any one species grows its biomass to a greater extent from externally supplied nutrients than from the excretions of other species. It does not mean, however, that cross‐feeding is unimportant for the assembly of species‐rich communities. This becomes clear from random viable metabolisms, whose environments provide only a single carbon‐source niche in the form of glucose. Without cross‐feeding, competitive exclusion dictates that only single species communities are possible in this environment, because a superior competitor would always replace the resident species (MacArthur & Levins, 1964; Stewart & Levin, 1973). In contrast, I find that communities with more than four species are formed even in this simple environment.

My predictions are consistent with recent experiments that cultured microbial communities from various natural habitats in a synthetic environment with only a single carbon source. Due to extensive cross‐feeding on metabolic by‐products, this simple environment could support communities of up to 22 species (Estrela et al., 2021; Goldford et al., 2018). I emphasize that this prevalence of competition over facilitation may depend on the environment and on the invading species. For example, in nutrient‐poor environments, or in environments where many species can grow only on metabolic excretions rather than externally supplied nutrients, facilitation may be more important than competition (Marsland et al., 2019).

Destructive invasions become more frequent during community assembly, but they never account for more than 50% of all invasions (Figure S2). This observation underlines that the external environment and metabolic excretions create ample free niche space to accommodate invaders without leading to the extinction of resident species. Constructive invasions, however, generally lead to smaller increases in community productivity than destructive invasions.

Genome‐scale models can provide a new perspective on long‐standing debates in invasion biology. One of these debates revolves around Elton's biotic resistance hypothesis. Evidence for this hypothesis is mixed, possibly because diversity may prevent invasions at small spatial scales, while facilitating invasions at large spatial scales (Altieri et al., 2010; Fridley et al., 2007; Jeschke et al., 2018). The present study system can help validate the hypothesis without space as a complicating factor, because it models a well‐mixed environment. And it shows that more species‐rich communities tend to be invasion resistant, because their efficient resource use leaves fewer opportunities for a new invader.

Elton argued that high species‐richness entails high invasion resistance. Part of the reason is that species‐rich communities may be associated with few ecological niches that are free for an invader to exploit. However, this association need not be universal. Many organisms create new ecological niches through “niche construction” or “ecosystem engineering” (Jones et al., 1994; Odling‐Smee et al., 2003). Bacteria do so by excreting by‐product metabolites that other bacteria can exploit through cross‐feeding. Such niche construction can help sustain large microbial communities even in chemically minimal environments (Estrela et al., 2021; Goldford et al., 2018; San Roman & Wagner, 2021). Thus, species‐rich communities might harbour more free niches, and these niches might also facilitate species invasions, contradicting Elton's reasoning. Indeed, my simulations show that species‐rich communities excrete a greater number of metabolites into the environment (Figure 3). The reason why they are nonetheless more invasion resistant is that competition for externally supplied nutrients prevails over facilitation via niche construction.

In the community assembly network I study, nodes are equilibrium communities. Two communities A and B are connected by a directed edge if a successful invasion of community A creates community B. Subsets of 100 species can form many more communities (≈1030) than one can explore. Even if only a small fraction of them form equilibrium communities in which all species coexist, the network of Figure 4a is only a sparse sample of the complete community assembly network. In addition, it is biased towards small communities, because I assemble successively larger communities from single species. Despite these limitations, the network representation can provide a new perspective on community assembly dynamics. For example, it leads to a natural distinction between “transient” communities, which have many outgoing edges and often give rise to other communities, and to “attractor” communities, which have many incoming edges and are frequent products of community assembly. Transient communities occur early during an assembly sequence, have low productivity and little invasion resistance. In contrast, attractor communities are late arising, with high productivity and high invasion resistance (Figure S9).

The most extreme attractor communities are those that have no outgoing edges – they are end points of assembly, because they cannot be invaded by any species. However, because the network of Figure 4a is a sparse sample of a complete assembly network, I cannot conclude that its 244 communities with outdegree zero are not invasible. This illustrates the limitations of sampling an assembly network. It could be overcome by studying assembly sequences for a much smaller species pool with a manageable number of possible communities. Complete assembly networks also could help answer a number of other questions that sparse sampling cannot. For example, they could help identify stationary distributions of community assembly, quantify non‐transitive interactions between communities and invading species through cycles in the network, and determine whether subsets of communities preferentially give rise to each other, thus creating assembly networks with a “modular” structure. Questions like these remain exciting tasks for future studies.

A main technical limitation of my work is its computational cost. For 250 assembly sequences of 5000 h (50,000 time steps of 0.1 h), maximal biomass‐growth needs to be computed 250 × 50,000 = 1.25 × 107 times, for each of a community's species, which is represented by a model of approximately 103 chemical reactions. Despite the numerical efficiency of FBA, this cost prevented me from asking a number of additional questions. For example, how does invasion resistance relate to other measures of community resilience, such as the ability to survive perturbations in species biomass or environmental nutrients. Do communities show multiple equilibria and priority effects? These questions too remain tasks for future studies.

While genome‐scale models are an important step towards biologically realistic models of community assembly, they also have biological limitations. First, they do not represent interference competition, in which bacteria release toxins or antibiotics to slow the growth of competitors. Second, they are not suited to represent viral parasitism or predator–prey interactions, which do occur in bacterial communities (Pérez et al., 2016), Third, they neglect that biomass growth is not only influenced by nutrient uptake. For example, different microbial species achieve maximal growth at different pH values, and metabolic by‐products often change the pH of the environment (Ratzke et al., 2020; Ratzke & Gore, 2018). Fourth, they neglect intercellular communication mechanisms such as quorum sensing, which can help a community coordinate the exploitation of nutrients (Goo et al., 2015).

The ability to persist over time is advantageous for a community, because it allows a community to “outlive” other communities. On the level of individual organisms, advantageous properties usually arise through Darwinian evolution. However, Darwinian evolution requires selection among competing individuals in a population. Because no analogue of a population exists during the assembly of a single community, an alternative mechanism is needed to explain increased resistance to invasions or other perturbations (Wagner, 2005). This mechanism has also been called “sequential selection”, “selection through survival alone” or “systemic selection” (Borrelli et al., 2015; Doolittle, 2014; Lenton et al., 2018). In this form of selection, ecological systems which experience changes that increase their ability to persist become more prevalent over time, because they outlive other, more ephemeral systems (Wagner, 2005). In my study system, these events are invasions of highly productive species, which increase the community's ability to resist future invaders. The same principle has been invoked to explain mechanisms that reduce oscillations in predator–prey systems, which make them prone to species extinctions (Doolittle, 2014), as well as planetary‐scale homeostatic mechanisms that are required for the Gaia hypothesis (Doolittle, 2014; Lenton et al., 2018). Darwinian evolution is not the only process that can create living systems with advantageous traits. The same can be achieved by ecological processes as fundamental as resource competition and repeated invasion.

AUTHOR CONTRIBUTIONS

The sole author performed all aspects of this work.

CONFLICT OF INTEREST

He declares that he has no conflict of interest pertaining to this work.

Appendix S1

Click here for additional data file.

Table S1

Click here for additional data file.