Neutrality in the Metaorganism

Michael Sieber; Lucía Pita; Nancy Weiland-Bräuer; Philipp Dirksen; Jun Wang; Benedikt Mortzfeld; Sören Franzenburg; Ruth A. Schmitz; John F. Baines; Sebastian Fraune; Ute Hentschel; Hinrich Schulenburg; Thomas C. G. Bosch; Arne Traulsen

doi:10.1371/journal.pbio.3000298

Abstract

Almost all animals and plants are inhabited by diverse communities of microorganisms, the microbiota, thereby forming an integrated entity, the metaorganism. Natural selection should favor hosts that shape the community composition of these microbes to promote a beneficial host-microbe symbiosis. Indeed, animal hosts often pose selective environments, which only a subset of the environmentally available microbes are able to colonize. How these microbes assemble after colonization to form the complex microbiota is less clear. Neutral models are based on the assumption that the alternatives in microbiota community composition are selectively equivalent and thus entirely shaped by random population dynamics and dispersal. Here, we use the neutral model as a null hypothesis to assess microbiata composition in host organisms, which does not rely on invoking any adaptive processes underlying microbial community assembly. We show that the overall microbiota community structure from a wide range of host organisms, in particular including previously understudied invertebrates, is in many cases consistent with neutral expectations. Our approach allows to identify individual microbes that are deviating from the neutral expectation and are therefore interesting candidates for further study. Moreover, using simulated communities, we demonstrate that transient community states may play a role in the deviations from the neutral expectation. Our findings highlight that the consideration of neutral processes and temporal changes in community composition are critical for an in-depth understanding of microbiota-host interactions.

Citation: Sieber M, Pita L, Weiland-Bräuer N, Dirksen P, Wang J, Mortzfeld B, et al. (2019) Neutrality in the Metaorganism. PLoS Biol 17(6): e3000298. https://doi.org/10.1371/journal.pbio.3000298

Academic Editor: Jeff Gore, MIT, UNITED STATES

Received: August 9, 2018; Accepted: May 15, 2019; Published: June 19, 2019

Copyright: © 2019 Sieber et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All relevant data are within the paper and its Supporting Information files (S1–S10 Data). Source code is available at https://github.com/misieber/neufit.

Funding: This research was supported by the Collaborative Research Centre 1182, "Origin and Function of Metaorganisms," funded by the Deutsche Forschungsgemeinschaft (https://www.dfg.de). TCGB gratefully appreciates support from the Canadian Institute for Advanced Research. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Abbreviations: AIC, Aikake information criterion; AMP, antimicrobial peptide; OTU, operational taxonomic unit

The microbial communities living in and on animals can affect many important host functions, including metabolism [1–3], the immune system [4–6], and even behavior [7, 8]. The extent and direction of this microbe-mediated influence is often linked to the presence or absence of species and their relative abundances. It is thus paramount to understand how host-associated microbial communities are assembled. A classical explanation for the emergence of a particular ecological community structure posits that every species is defined by distinct traits and occupies a specific ecological niche. An implicit hypothesis underlying much of microbiome research is that hosts have the potential to actively shape their associated microbial communities by providing niches for useful microbes [9]. This implies that individual hosts could select for a potentially very specific community structure.

While the assumption that the metaorganism—the host together with its associated microbes [10, 11]—is an actively shaped symbiotic unit is appealing, it may bias interpretation of microbiota-host analyses. An example is the widely reported connection between the structure of the human gut microbiota and obesity, which turns out to be very hard to distinguish from random noise [12]. More recently, it has also been demonstrated that the genetic background of the host does not significantly shape human microbiome composition [13]. Indeed, it has been argued that selective processes should not form the null hypothesis for explaining the allegedly cooperative host-microbe symbiosis [14]. Instead, neutral models have been proposed as a valuable tool for finding patterns in the tremendous complexity of ecological communities that may have a deeper mechanistic cause. Neutral models assume ecological equivalence between species. Thus, the community structure within a single host is the outcome of purely stochastic population dynamics, immigration, and local extinctions [15]. The Unified Neutral Theory of Biodiversity [16, 17] extends the neutral framework from a single site to multiple sites, each harboring its own local community (Fig 1). Diversity within the local communities is maintained by immigration from a common source community, so that the whole setup resembles a mainland-island structure [18]. Neutral theory has been applied to numerous ecological systems and sparked a lot of controversy along the way [19].

Download:

Fig 1. Sketch of the setup of the neutral model.

Only a subset of the available microorganisms can pass through the environmental filter, forming a common source pool for all hosts. This selective step is not described by the neutral model. The microbiota of individual hosts form local communities, and the microbial intra-host population dynamics are neutral. Local diversity is maintained by immigration from the source pool.

https://doi.org/10.1371/journal.pbio.3000298.g001

Despite these controversies and with the growing availability of microbial community data, the neutral theory has been adapted and applied to the microbial world [20–24]. In particular, the local community versus metacommunity structure of neutral theory naturally extends to host-associated microbial communities. Here, the hosts are viewed as ecosystems and their microbiota are treated as local communities [11, 25]. This has led to several recent studies addressing the question of whether the microbiota conform to the patterns predicted by the neutral theory. The rank-abundance patterns of the microbiota from three domesticated vertebrates, for example, were found to be largely consistent with the neutral expectation [26]. A good fit of the neutral model was also found for the microbiota of young individuals of the zebrafish Danio rerio [27]. Interestingly, in this study neutrality decreased as hosts aged, indicating an increasing influence of selective processes over developmental time, potentially linked to the activation of a fully functioning adaptive immune response. The human microbiota were found to be predominantly non-neutral across various body sites [28], which is thought to be the result of a phase transition from a dispersal-dominated neutral regime to a selection-dominated within-host regime [29]. In this study the few neutral communities were mostly associated with the urological tract and the skin, while in another study the composition of the skin microbiota of healthy human subjects in large Chinese cities was found to be better explained by non-neutral processes [30]. Contrasting results can also be obtained depending on the state of the host, with the healthy human lung microbiota being largely consistent with a neutral model, while microbes recovered from diseased lungs diverged from neutrality [31]. Neutral assembly processes in invertebrate hosts are understudied, but a recent study shows the microbial communities associated with the fruit fly Drosophila melanogaster to be consistent with the predictions of a neutral model [32]. While these studies have drawn awareness to the potentially important contribution of neutral processes in shaping the microbiota, the use of several different neutral models and focus on one or a few, mostly vertebrate host organisms, makes it hard to draw more general conclusions.

Here, we aim to overcome these limitations by consistently applying a neutral model to a variety of distinct host systems, including but not limited to a wide range of invertebrate hosts. We included host species from a range of eukaryotic multicellular organisms with different lifestyles, ranging from early branching groups such as sponges to the house mouse, with its fully developed adaptive immune system. We used a general modelling approach and combined it with a consistent simulation framework in order to achieve an unbiased comparison of the distinct metaorganisms. Based on the neutral null hypothesis, we identified members of the microbiota that consistently deviate from the neutral expectation, possibly indicating differential selection of these particular microbes. Additionally, we hypothesize that some instances of observed non-neutrality may be due to transient community stages, in which the microbiota has not yet reached its long-term equilibrium composition.

Results

We employed the neutral model presented by Sloan and colleagues [20], which is particularly suited for large-sized microbial populations and has been used before to assess microbiota neutrality [27, 32]. It describes the stochastic population dynamics within a local community, corresponding to the microbiota of a specific host. To maintain local diversity, which would otherwise reduce to a single species through ecological drift, immigration from a fixed source community occurs with rate m. This source community is not equivalent to the environmental pool of microorganisms, but rather it can be interpreted as the collection of all microbial species that can pass through the environmental filter of the host (Fig 1). Thus, the neutral model does not make any assumptions about the selective constraints that apply before and during colonization of the host and, in particular, it is not concerned with any host traits that may restrict the range of colonizing microbes [33].

This model allows to derive an expression for the expected long-term stationary community composition. Specifically, it predicts the relationship between the mean relative abundance of a taxon across all local communities, i.e., the metacommunity, and the probability of actually observing this taxon in any single community. It has been shown that this relationship is determined by a beta distribution ([20] and Methods). The only free parameter of this model is the immigration rate, m, which can be calibrated by a nonlinear fit. The goodness of the fit then indicates how well the prediction of the neutral model compares to the empirical data. More specifically, we used the coefficient of determination R² as a quantitative measure of how consistent the data are with the neutral model. See Methods for details of the neutral model, simulations, and the fitting procedure.

Neutrality of the microbiota

We fitted the theoretical neutral expectation to the published microbiota compositions of eight different host species across four phyla and, for comparison, three environmental microbial communities (see Methods and S1 Table for an overview and references).

With the species Sarcotragus fasciculatus, Ircinia oros, and Carteriospongia foliascens, our study included examples from the oldest extant sister group to all other animals, the sponges (Porifera) [34]. Despite their filter-feeding lifestyle, sponges harbor a very diverse and highly specific microbiota [35], which mediates the functional role of the sponges in the ecosystem [36]. This is complemented by the jellyfish Aurelia aurita, the starlet sea anemone Nematostella vectensis, and the freshwater polyp Hydra vulgaris. These hosts are examples from another early branching phylum, the Cnidaria, which possesses a basic innate immune system thought to play an important role in controlling the interaction with the associated microbial community [37]. We also included the nematode Caenorhabditis elegans as one of the best-studied multicellular organisms, whose microbiota has only more recently come into the focus [38–41]. Finally, the house mouse Mus musculus with its potent adaptive immune system [42] offers a well-studied intestinal microbiota [43, 44]. Our datasets included lab-reared animals as well as samples from several natural populations from various locations, representing a broad range of environmental and microbial contexts. For some host species, samples from both lab and natural populations were available, which we analyzed separately to assess the potential impact on neutrality of heterogeneous natural habitats versus homogeneous, constant lab conditions. See Methods for a summary of the context of the individual datasets.

In all cases, community composition was determined by the relative abundances of operational taxonomic units (OTUs), obtained by standard high-throughput 16S rRNA sequencing techniques [45]. To ensure equal sample sizes, all OTU abundance tables were rarefied to the same read depth (1,000 reads per sample). We analyzed the effect of the rarefaction on the neutral fit for datasets in which more reads were available and found it to be of little relevance for most communities (S6 Fig).

A comparison of the theoretical and observed relationship between mean relative abundances and occurrence frequencies of OTUs for three illustrative examples is shown in Fig 2 (see S1 Fig for examples of the neutral fit for all datasets). The neutral model generally predicts a specific monotonic increase of the occurence frequency of an OTU, with an increasing mean relative abundance of this OTU across all samples. This is a quantitative reflection of the null expectation that more abundant OTUs should also be found in more samples. Consequently, OTUs that are found below the prediction in the lower right region of the graph are found in fewer samples than expected by their mean abundance across all samples. Conversely, OTUs that lie above the neutral expectation are found more often than expected.

Download:

Fig 2. The relationship between mean relative abundance of taxa across all samples and the frequency with which they are detected in individual samples for three representative datasets.

Each dot represents a taxon and the solid line is the best-fitting neutral community expectation. The dashed lines and gray area depict the 95% confidence bands. The left panel shows a microbial community sampled from a compost as an example of an abiotic environment. The center panel shows data for the natural microbiota of C. elegans nematode populations sampled from the same compost. The right panel shows the microbiota of a wild M. musculus mouse population. Here, the red dots indicate members of the microbial family Ruminococcaceae, which are predominantly overrepresented compared with the neutral expectation. A genus that was found in much fewer mice than expected by chance is Shigella, being present in less than a quarter of the hosts despite a neutral expectation of nearly 100% prevalence. The OTU abundance data used to generate this figure are available in S6 and S7 Data. OTU, operational taxonomic unit.

https://doi.org/10.1371/journal.pbio.3000298.g002

A summary of the goodness of fit of the neutral model to the microbial communities from all datasets is shown in Fig 3. Higher R² values generally indicate a closer match between the neutral expectation and the data. To aid the interpretation and comparison across datasets, we additionally contrasted every dataset with a corresponding simulated, completely neutral community that had the same sample size, microbial diversity, and estimated immigration parameter, m (see Methods for details of the simulation). This allows a comparison of the deviation between the empirically obtained measure of neutrality and the expected values from a known neutral community for each dataset. This is crucial because, depending on the sample size and microbial diversity, even a completely neutral community will not yield a perfect fit, even when there is sufficient time to reach an equilibrium (Fig 4).

Download:

Fig 3. Goodness of fit of the neutral expectation to a range of host-associated and environmental communities.

Circles denote natural populations and diamonds denote laboratory populations; error bars indicate 95% bootstrap confidence intervals. The gray bars indicate the 95% bootstrap confidence intervals around the expected goodness of fit (horizontal black line) to a neutral simulation, providing a comparison to a completely neutral community with the same sample size and diversity. The data for C. foliascens are from several different natural populations, while the data for H. vulgaris and N. vectensis are from different time points. The spread of points along the x-axis is added to increase visibility. The data used to generate this figure are available in S7 Data. Phylogeny were generated with phyloT based on NCBI taxonomy. Sponge photographs courtesy of Susanna López-Legentil (UNC Wilmington) and Mari-Carmen Pineda (AIMS). AIMS, Australian Institute of Marine Science; NCBI, National Center for Biotechnology Information; UNC, University of North Carolina.

https://doi.org/10.1371/journal.pbio.3000298.g003

Download:

Fig 4. Impact of sample number on the goodness of fit of the neutral expectation.

For each sampling depth, communities were resampled 100 times with replacement at a rarefaction level of 1,000 reads/sample. The left panel shows the results for simulated neutral communities. The immigration parameter was fixed at m = 0.05, and neutral dynamics were simulated for 10⁶ time steps in all cases. Generally, goodness of fit remained high over a large range of sample sizes, reflecting the underlying completely neutral community dynamics. Note that the goodness of fit did not reach 1.0 in finite and relatively small metacommunities, even if all communities were sampled. The right panel shows the effect of subsampling for several host-associated microbial communities. As for the simulated data, observed levels of neutrality were relatively robust, only dropping off at very low sample sizes. The C. elegans datasets showed a very poor fit at very low sampling depths but began to level off well below the neutrality values of the other datasets. Shaded areas indicate 95% bootstrap confidence intervals. The data used to generate this figure are available in S1 Data. OTU, operational taxonomic unit.

https://doi.org/10.1371/journal.pbio.3000298.g004

Turning first to the environmental microbial communities, we generally observed a high consistency with the neutral expectation. The majority of microbial taxa found in compost samples, for example, closely followed the neutral prediction (Fig 2), as reflected by a high goodness of fit measure (R² = 0.86). While the goodnesses of fit to the microbial communities in seawater and marine sediment are lower, they lie in the range expected from the respective neutral simulations. The sediment dataset showed a marked effect of rarefaction, with the goodness of fit increasing to the same level as the seawater community when more reads are included (S6 Fig). It is also the dataset with the highest microbial diversity (>3,000 OTUs) and a relatively low sample size, which may indicate a greater influence of rarefaction on rare taxa. The overall good consistency of the environmental communities with the neutral expectation aligns with previous studies showing that random population dynamics and immigration play a major role in shaping microbial communities in abiotic environments [21, 23, 24].

Biological hosts, with their potential for active selection of specific taxa, could on the other hand be expected to have their resident microbial communities under tighter control, reducing the relative importance of stochastic processes after colonization, which would be reflected in a worse fit of the neutral expectation. The nematode C. elegans, for example, naturally consumes bacteria as its food source and possesses an innate immune system as a defence against the threat of ingested pathogenic microbes [46]. This suggests that these worms have some control over their microbiota, and indeed we found a substantial mismatch between the neutral expectation and the worm microbiota (Fig 2), resulting in a low goodness of fit (R² = 0.44). We also found a similar level of neutrality for the microbiota of laboratory worm populations (Fig 3), indicating a similar influence of neutral processes under more controlled laboratory conditions. For both the natural and lab-enriched worms, the goodness of fit was well below the levels expected from the corresponding neutrally simulated communities (Fig 3). It is noteworthy that these specific worms were isolated from the same compost where the environmental microbial communities showed a very good alignment with the neutral expectation (Figs 2 and 3).

Vertebrate hosts, with their more sophisticated and complex immune systems, are expected to diverge even further from the neutral expectation. However, intriguingly, the neutral model showed a very good fit to the microbiota of a wild M. musculus mouse population (Fig 2). The best fit was in fact on par with the fit to the environmental microbes isolated from compost (R² = 0.86), and the deviation from the simulated neutral community was minimal (Fig 3). The high level of neutrality in the natural mouse population was corroborated by a laboratory population of M. musculus (Fig 3).

We found similar high levels of neutrality for sponges, which, as aquatic filter feeders, are exposed to the full range of marine microbes. Like nematodes, they have a basic immune system [47], and, especially for S. fasciculatus and some populations of C. foliascens, the data were in very good agreement with the neutral expectation (Fig 3 and S1 Fig). In contrast, the results for the three aquatic polyps were less consistent. The two laboratory populations of A. aurita were in very good agreement with the neutral prediction, while samples taken at different developmental stages of H. vulgaris and N. vectensis generally showed a larger deviation from full neutrality (Fig 3).

The estimated immigration parameter, m, which can be interpreted as a proxy for inter-host transmission and similarity of the individual microbial communities, showed considerable variation across the datasets (S2 Fig). Overall, communities from aquatic samples showed higher best-fit values of m than terrestrial samples. This potentially indicates an increased between-sample dispersal and a more homogeneous microbial distribution in aquatic settings, at least on the relevant spatial scales. This seems especially plausible for the filter-feeding sponges sampled from the same location. In contrast, we obtained relatively low estimates for the immigration parameter for the mice and nematodes, in both the natural and the laboratory populations. For nematodes, in particular, a more patchy distribution of microbes on the spatial microscale, together with stochastic founder effects during colonization [48], may play a role in creating more heterogeneous microbiota compositions. In combination with less exchange of microbes in this terrestrial organism with low dispersal ability, this can contribute to low values of m.

Note that we did not find a correlation between immigration rate m, sample size or microbiota diversity, and the level of neutrality (S5 Fig). We additionally analyzed the effect of varying sample sizes by subsampling all datasets with more than 20 samples, because the number of individual samples has the potential to influence the composition of the source community of available microbial species. We found the observed levels of neutrality to be robust against subsampling, unless sample sizes are very small (Fig 4). In all cases, neutrality consistently decreased with decreasing sample sizes, potentially reflecting a decreasing overlap in the local microbial communities. We performed this analysis for both the simulated and the empirical datasets, which both showed a very similar impact of subsampling.

Comparing the neutral expectation with data from a wide range of environments and host species showed that a good fit of the neutral model is in fact the norm rather than the exception (Fig 3). In particular, there was no substantial difference in observed neutrality between microbial communities living in environmental samples, such as compost and seawater, and the microbiota from several animal host species. We thus conclude that the observed patterns of microbiota composition in hosts as different as mice and marine sponges can in many cases already be obtained with a neutral null model.

Identifying non-neutral taxa

Our analysis showed that for all host species, regardless of the specific goodness of fit of the neutral model, there are microbial taxa that diverged substantially from the neutral expectation. Because an essential ingredient of neutral theory is the assumption of selective equivalence, such a divergence from neutrality may signify taxa that are subject to differential (e.g., host-mediated) selection. Thus, identifying specific non-neutral members in the bulk of the microbiota is of great practical importance. In particular, OTUs that are observed more often than expected by chance are commonly seen as good candidates for beneficial symbiotic microbes. In general, however, it is not straightforward to decide whether a microbial taxon is under positive or negative selection within a specific host. A taxon could, for example, be observed more often than expected because it is indeed essential for the host, even at low abundances, and is thus positively selected for in the community. Alternatively, such overrepresentation may point to a microbe that is simply a very good colonizer, i.e., it can easily disperse and pass through a host's environmental filter, but within the host it is subject to negative selection to prevent it from reaching high abundances. Here, the practical value of the neutral model lies in its ability to identify microbes that appear to interact in a specific, non-neutral form with the host and are thus worthy of further investigation, taking into consideration ecological and functional information.

In the following, we classified all OTUs outside of the 95% confidence bands around the neutral prediction as not consistent with the neutral expectation. We found that the distribution of non-neutral OTUs around the neutral prediction is roughly symmetrical in all cases, and no microbial order deviated systematically above or below the neutral expectation (S3 Fig).

To further distinguish random deviations from neutrality and actual non-neutral candidates, we also applied a more conservative definition of non-neutrality. For this, a particular OTU was classified as non-neutral only when it consistently diverged from the neutral expectation in the same direction in independent host populations. Applying this definition to the natural and laboratory populations of C. elegans and M. musculus yielded a reduced subset of microbial taxa, which lay above or below the neutral prediction in both populations.

This analysis revealed that, of the C. elegans microbiota, only the genus Ochrobactrum is consistently underrepresented (S2 Table). It is found in only about 40% of the natural isolates and 70% of the laboratory worms, despite its relatively high mean abundance across all worms and a neutral expectation of 100% prevalence. This underrepresentation of Ochrobactrum is intriguing, as it had previously been identified to be enriched in worms, to favor growth of worm populations, and to be able to persist in the nematode's gut even under starvation conditions [39]. However, this observation is consistent with a transient community state (see below), where Ochrobactrum has not yet colonized all worms, despite being very successful once it has entered the worm.

In the M. musculus microbiota of natural and laboratory populations, over a third of the overrepresented taxa were from the bacterial family Ruminococcaceae. Members of this family were never consistently found below the neutral expectation across hosts (Fig 2 and S3 Table). This family includes several physiologically relevant genera involved in the utilization of plant polysaccharides such as starch and cellulose [49], and high-fat diets low in plant-derived materials have been found to decrease the proportion of the Ruminococcaceae in mouse guts [50]. A notable microbe that was found in significantly fewer mice than expected from the neutral model is Shigella (Fig 2). Invasive Shigella can cause intestinal inflammation in humans, but mice mount an effective defense against it, thereby preventing acute infections [51]. This illustrates how potentially pathogenic microbes, which may be actively selected against by the host, can be found diverging from the neutral expectation. That we were able to consistently identify a bacterial family that had previously been linked to important metabolic functions as overrepresented and potential pathogens as underrepresented shows the utility of the neutral model as a null hypothesis and a tool of finding meaningful patterns in the vast complexity of the microbiota.

Another commonly applied approach to tackle this complexity consists of the identification of OTUs that are shared among hosts, typically denoted as the core microbiota. This usually requires setting a specific occurence threshold for the taxa found across hosts in similar habitats [52]. Applying this persistence-based definition to our data shows that a substantial portion of the core taxa, e.g., those with an occurence frequency above 90%, fall within the neutral expectation (Fig 2 and S1 Fig). This implies that many taxa within the abundant and prevalent core microbiota may actually not occur more often than would be expected by chance. Taxa such as the Ruminococcaceae in mice, on the other hand, which were below typical cutoff prevalences but occurred more often than expected by chance, emerge as interesting candidates beyond the usual core taxa. Deviations from the neutral expectation can thus complement and enrich insights gained from the widely used concept of the core microbiome.

Change of neutrality over time

As we have discussed in the previous section, the most commonly assumed cause of deviations from neutrality is selection by the host, for example, through its immune system. In the context of our study, this would mean that after colonization, host-mediated selection may have altered the composition of the microbial community, resulting in a non-neutral community structure. Interestingly, we found a particularly strong deviation from neutrality for the nematode host, which possesses a rather simple immune system [53]. In contrast, other hosts with more complex immune systems, for which strong selection could have been presumed, show microbiomes close to neutral expectations. This suggests a role of additional factors beyond the immune system driving deviations from the neutral expectation.

A crucial aspect of the neutral prediction is that it is based on the expected long-term stationary distribution of the microbial community. If the community is in a transient state, however, it can appear non-neutral even if it is the product of a purely neutral process. This effect is illustrated in Fig 5, showing a simulation of the population dynamics of 50 neutral communities (Methods). For every community, only a few random colonizers were picked from the source community, which yields random but non-neutral initial communities. We then compared the neutral expectation to the in silico communities over time during the simulation, showing that they become "neutralized" through random population dynamics and immigration until the community composition reaches its long-term equilibrium, or the host dies.

Download:

Fig 5. Goodness of fit of the neutral long-term expectation to neutrally simulated communities over time.

The two insets show the best fits of the neutral model to snapshots of the communities at an early time point (indicated by the red dot) and a late time point (indicated by the blue dot). Scale and meaning of the axes for the insets are the same as in Fig 2. Each community was initialized with a few random colonizers, which yielded non-neutral initial communities that neutralized over time through a stochastic death-birth process and immigration. The shaded area indicates the 95% bootstrap confidence intervals. The data used to generate this figure are available in file S1 Data.

https://doi.org/10.1371/journal.pbio.3000298.g005

Thus, if community assembly is governed by purely neutral processes, we would expect the fit of the neutral model generally to become better over time. Decreasing neutrality, which has been reported for the zebrafish, D. rerio [27], would on the other hand be a good indicator of selective processes. A similar increase in the relative importance of selective processes over stochastic effects has been described for the ecological succession of microbial communities in salt marshes [54]. However, while the microbiota from the Hydra and Nematostella populations show clear developmental patterns [55], we did not find a clear temporal trend in the deviation from the neutral model (S4 Fig).

While we could not conclusively disentangle transient effects in the interplay between ecological drift and selective processes, in principle, it is possible that a poor fit of the neutral expectation may be due to a transient non-equilibrium state of the community rather than actual non-neutral processes. This effect would be especially pronounced for short-lived host species such as C. elegans, in which initial colonization can lead to communities that appear non-neutral even for ecologically equivalent (i.e., neutral) colonizers [48]. The microbiota of such hosts may never effectively converge to the neutral long-term expectation within the lifetime of the host, even if they are dominated by stochastic dynamics and immigration. For longer-lived species, such as sponges and mice, we would expect these transient effects to play a smaller role.

Comparison with a random sampling model

The observed high levels of neutrality do not necessarily imply that the microbiota is simply a random collection of microbes. A good fit of the expected distribution indicates that the simple neutral model is sufficient to describe the relationship between abundance of a species and its occurrence frequency. But this cannot rule out other processes, niche related or stochastic, that lead to the same community structure [56].

A very simple alternative model is obtained by assembling a local community by drawing randomly from the fixed source community, thus neglecting stochastic reproduction events and dispersal. In this case, the probability of observing a specific species in a local community is determined by a binomial distribution. Comparing the best fits of the neutral expectation and of a binomial distribution using the Aikake information criterion (AIC) indicated that the neutral model fits the data better in all cases (S1 Table). This indicates that the microbiota are not merely random samples of the microbes present in the metacommunity.

In a more general context, ecological communities such as microbiota are examples of component systems, in which each specific instance of the system (e.g., a specific microbial community) is composed of components drawn from a shared set of basic building blocks (e.g., the source pool in our scenario). Such systems have been shown to exhibit some general invariant properties, and deviations from these properties may reveal functional constraints imposed on the shared components [57].

Discussion

We set out to quantify how consistent the microbiota of a range of different host organisms are with the expectation from a neutral model. We found that the structures of the considered host-associated microbial communities are often in surprisingly good agreement with the neutral expectation. Taking into account a more complete simulation analysis of the neutral model further suggested that transient stages of the microbiota may contribute to the divergence from neutrality observed for some host species.

Even a high consistency with the neutral expectation does not rule out the presence of selective abiotic or biotic processes, as non-neutral models can yield predictions consistent with neutral theory under certain conditions [56]. Moreover, biological hosts will select for a subset of the environmentally "available" microbes, and yet the resulting microbiota may still look mostly neutral. In fact, comparing the microbes present in seawater to the species present in individual sponges reveals that there is almost no overlap (S7 Fig). This indicates that sponges are a highly selective environment, imposing strong constraints on potential colonizers, and yet after colonization, their microbiota are very consistent with the neutral expectation. A similar selective filter is found for C. elegans, in which only a fraction of the microbes present in the environment successfully colonize the worms (S7 Fig). This suggests that there are host traits that act as environmental filters [33], and they often will do so in a predictable, i.e., deterministic way. In this sense, the microbiological tenet that "everything is everywhere, but, the environment selects" [58] also applies when the environment happens to be another biological organism. However, our results highlight that the processes that are at play after the environmental filter has been passed can be highly consistent with the neutral model. This apparent contradiction between selective environments and neutral processes acting on the allowed set of species may also have contributed to some of the contrasting results found in previous studies [59].

This filtering of potential colonizers is not unique to biological hosts, but what sets biological hosts apart from other habitats is the potential of the hosts' filtering mechanisms to evolve. An example for such an evolvable mechanism is the production of antimicrobial peptides (AMPs) by several Hydra species [60]. Here, each Hydra species was shown to possess a specific AMP expression profile, selecting for a species-specific set of bacteria by inhibiting colonization of foreign microbes. Another recent study showed how host-derived proteins from the jellyfish A. aurita interfere with bacterial quorum sensing to modulate host colonization [61]. While these examples do not constitute proof of adaptive coevolution between the host and its microbiota, they exemplify the possible specific selection capabilities of the host. Moreover, our results then highlight that microbiota composition is the outcome of several intertwined processes. For example, selective filtering at the colonization stage may be controlled by the host or by microbe characteristics, and is then presumably followed by largely neutral, stochastic processes. It is further possible that selection at either the colonization stage or community assembly thereafter only affects a subset of microbes, while the remaining microbes follow neutral dynamics. Another factor to consider is that taxonomic and functional microbial community composition may be shaped by largely independent processes [62, 63]. Given the known examples of functional redundancy present, e.g., in the human gut microbiota [64], viewing communities from the perspective of functional (e.g., metabolic) categories can provide an alternative picture to data based on 16S rRNA gene-based taxonomic profiling [65]. These intertwined processes need to be better understood and taken into account for an enhanced understanding of microbiota assembly.

Our simulation analysis demonstrates that the community needs time to reach the dynamic neutral equilibrium (Fig 5) through ecological drift. This emphasizes the need for temporal data to further investigate the role of transient community stages and drift in microbial communities. Such data are usually not available for host-associated microbial communities, and our study highlights the need to be aware of this potentially confounding factor when applying neutral models to cross-sectional community snapshots.

The neutral assumption that species do not interact is anathema to many ecologists, and for the microbiota in particular it is assumed that the functioning of the community is enhanced by cooperative interactions. Recent results, however, suggest that interactions within the mouse gut microbiota are indeed predominantly competitive and very weak [66]. Together with our finding of the microbiota often being consistent with neutral expectations for several different host species, this lends support to taking the neutral null hypothesis as a key component of our understanding of host-microbe interactions. We would like to emphasize that even a good consistency with the neutral model does not imply that the microbiota is functionless or that hosts are nonselective—the Hydra system actually contains taxa that confer a selective advantage, but which appear as neutral in our present analysis [67]. Neutral theory in fact makes no claims about microbiota function; rather, it assumes that the different community compositions that arise after the environmental filter of the host has been passed are, in the words of the original formulation of the neutral theory, selectively nearly equivalent, that is, they can do the job equally well in terms of survival and reproduction [68] of the individual host.

Awareness of the potentially significant role of neutrality in shaping the microbiota is crucial to avoid being led astray by randomness in view of the incredible complexity of host-microbiota symbioses and to be able to uncover general principles of microbiome assembly.

Methods

Summary of datasets and data availability

The collection of microbial abundance datasets for the environmental and host-associated communities (i.e., OTU tables) was assembled as part of the Collaborative Research Centre 1182 "Origin and Function of Metaorganisms" (https://www.metaorganism-research.com) funded by the Deutsche Forschungsgemeinschaft. All OTU tables are available in the Supporting information accompanying this paper (S2–S10 Data). The raw sequencing data have been published previously and are available online; see below and S1 Table for the corresponding references. In the following, we give a short summary of the context of the respective datasets; for details we refer to the corresponding references.

Sponges [35].

S. fasciculatus (S2 Data) and I. oros (S3 Data) samples were collected at the same location on the Spanish coast on different dates. C. foliascens samples were collected at different locations along the Australian coast on different dates (S4 Data).

C. elegans [39].

The natural isolates of single C. elegans worms were collected on different dates at the same location in Northern Germany (S6 Data). For the lab-enriched populations of C. elegans (S6 Data), single worms isolated from the same location were transferred to agar plates seeded with Escherichia coli as food. On these plates, the worms were allowed to reproduce via selfing and were maintained for 2–3 weeks before sampling the resulting populations.

Mice [44].

For the natural mice, 69 mice were captured from 34 unique locations, with a maximum of 3 mice per location, in the hybrid zone of M. musculus musculus and M. musculus domesticus in Bavaria, Germany, during May and June 2011 (S7 Data).

The lab dataset includes 40 F2 hybrids between M. musculus domesticus and M. musculus musculus standard lab strains and, additionally, 7 each of the M. musculus domesticus and M. musculus musculus parental mice (S7 Data). The parental strains were kept under conventional conditions for two generations before setting up experimental crosses.

N. vectensis [55].

Animals were sampled from 15 independent laboratory cultures under different but fixed environmental conditions (temperature and salinity) at specific time points (S8 Data).

H. vulgaris [69].

Animals were sampled at several time points from independent, clonal lab cultures with identical and constant environmental conditions (S9 Data).

A. aurita [61].

Animals were sampled from 4 independent populations: 1 control population and 3 treatment populations (S10 Data). The treatment (QQ) populations were each treated with a different host-derived quorum-quenching protein.

Environmental samples.

Compost samples were collected in the context of the C. elegans study [39], following the same protocols and at the same locations (S6 Data). Marine seawater and sediment samples were collected in the context of the sponge microbiota study [35] (S5 Data).

The neutral model

Here, we summarize the model presented in [20], considering only the purely neutral special case with equivalent growth rates for all species. A local community contains a fixed number of N individuals, of which N_i are of species i. Individuals die with a constant death rate δ, and a dead individual is replaced by reproduction of a local individual (with probability 1−m) or by immigration from a fixed source community (with probability m). Thus, the total population size N remains constant. The relative abundance of a species in the source community is denoted by p_i. With this, the probabilities that the number of individuals of species i increases by one, decreases by one, or stays the same are given by (1) These transition probabilities correspond to Hubbell's original neutral model, but the following continuous approximation derived by Sloan and colleagues [20] can be efficiently applied to very large population sizes and allows for a relatively simple analytical solution. In particular, this allows for a calibration with the high-throughput 16S rRNA sequencing data obtained from microbial communities.

If N is large enough, such as in microbial communities, we can assume the relative abundance x_i = N_i/N of species i to be continuous. This leads to a Fokker-Planck equation for the probability density function ϕ_i(x_i,t) of x_i, (2) where M_δx and V_δx are the expected rates of change in relative abundance and variability, respectively. Assuming the time intervals between individual death-birth events are short, they can be approximated by (3) and (4)

Eq 2 describes the neutral population dynamics of a local microbial community. This is not, in general, amenable to analytical treatment, but one can approximate the long-term equilibrium solution of this equation. Namely, the so-called potential solution for ∂ϕ_i/∂t = 0 arising from is given by the beta distribution (5) with c = Γ(Nm)/[Γ(Nm(1−p_i))Γ(Nmp_i)]m. This can be connected to empirical observations. For a given detection threshold d, the probability of actually observing a species in a local community is given by the truncated cumulative probability density function (6)

Here, N is given by the number of reads per sample, and the relative abundance p_i of the focal species in the source community can be approximated by the mean relative abundance of the species across all samples. The detection threshold is set to d = 1/N. The probability of immigration m is thus the only free parameter and can be used to fit the predicted long-term distribution to the observed occurence frequencies of the microbiota.

Fitting procedure

The fitting process applied to all datasets, empirical and simulated, works as follows. First, the OTU abundance table was rarefied to the desired read depth (N = 1,000 reads per sample for the main results, which was determined by the lowest available read depth from all datasets). The expected observation probability (Eq 6) and a binomial model were then fitted to the observed mean relative OTU abundances and occurence frequencies f_i obtained from this rarefied table using nonlinear least squares minimization with the lmfit package (pypi.org/project/lmfit) [70]. As a measure of the goodness of fit, we then calculated the standard coefficient of determination using the ratio of the sum of squared residuals and the total sum of squares: (7) where Φ_i is the expected ocurrence frequency obtained from the best-fit neutral prediction. Additionally, for comparison of the neutral and binomial model, the mean AIC was calculated. In all cases, 95% bootstrap confidence intervals were obtained by resampling the hosts 100 times with replacement and performing the described fitting procedure on each resampled set of hosts.

Neutral model simulations

Simulated neutral communities were generated following the stochastic death-birth-immigration process described above and in [20]. Communities with a prescribed number of samples/communities, community size, source pool diversity, and immigration parameter m were started from random initial conditions and simulated for 10⁶ time steps to let them reach dynamic equilibrium. The community snapshots obtained at the end of each simulation were treated with the same fitting procedure as the empirical datasets.

Code availibility

The Python code for fitting and simulating the neutral model is available at https://github.com/misieber/neufit.

Supporting information

S1 Table. Overview of the datasets included in this study, their source (natural or laboratory communities), number of samples, microbial diversity, and a summary of the results of the model fits.

Best-fit data are the mean values from 100 rarefactions at a level of 1,000 reads per sample. Reference is to the original study containing detailed information on sample preparation and raw sequencing data. The data summarized in this table are also available in S1 Data.

https://doi.org/10.1371/journal.pbio.3000298.s001

(PDF)

S2 Table. Non-neutral genera found in the C. elegans microbiota.

https://doi.org/10.1371/journal.pbio.3000298.s002

(PDF)

S3 Table. Non-neutral genera found in the M. musculus microbiota.

https://doi.org/10.1371/journal.pbio.3000298.s003

(PDF)

S1 Fig. The relationship between mean relative abundance of taxa across all samples and the frequency with which they were detected in individual samples.

Each dot represents a taxon and the solid line is the best-fitting neutral community expectation. The dashed lines and gray area depict the 95% confidence bands. The OTU abundance data used to generate this figure are available in S2–S10 Data. OTU, operational taxonomic unit.

https://doi.org/10.1371/journal.pbio.3000298.s004

(PDF)

S2 Fig. Overview of the best-fit values for the immigration parameter m of the neutral model.

Circles denote natural populations and diamonds denote laboratory populations; error bars indicate 95% bootstrap confidence intervals. The data for C. foliascens are from several different natural populations, while the data for H. vulgaris and N. vectensis are from different time points. Spread of points along the x-axis is added to increase visibility. The data used to generate this figure are available in S1 Data. Phylogeny were generated with phyloT based on NCBI taxonomy. Sponge photographs courtesy of Susanna López-Legentil (UNC Wilmington) and Mari-Carmen Pineda (AIMS). AIMS, Australian Institute of Marine Science; NCBI, National Center for Biotechnology Information; UNC, University of North Carolina.

https://doi.org/10.1371/journal.pbio.3000298.s005

(PDF)

S3 Fig. Bars show the mean relative abundance of the 10 most abundant microbial orders in a specific host population.

The colored sections indicate the fractions of OTUs within that order that were found above the neutral expectation (red), within the neutral expectation (gray), and below the neutral expectation (blue). The OTU abundance data used to generate this figure are available in S2–S10 Data. OTU, operational taxonomic unit.

https://doi.org/10.1371/journal.pbio.3000298.s006

(PDF)

S4 Fig. Neutrality of the microbiota of H. vulgaris and N. vectensis at different time points during development and adult life stages (95% bootstrap confidence intervals attached).

The dashed lines are the best linear fits; slopes are not significantly different from zero. The data used to generate this figure are available in S1 Data.

https://doi.org/10.1371/journal.pbio.3000298.s007

(PDF)

S5 Fig.

Consistency with the neutral model versus the estimated dispersal parameter m (top left), the number of samples (top right), the number of identified taxa (bottom left), and the Shannon index of diversity (bottom right). Circles denote natural populations and diamonds laboratory populations. The data used to generate this figure are available in S1 Data.

https://doi.org/10.1371/journal.pbio.3000298.s008

(PDF)

S6 Fig. Impact of rarefaction on neutrality.

Generally, consistency with the neutral model initially decreased with increasing read depth, until it leveled off as read depth increased further. For communities that showed a high consistency with the neutral model, varying read depths did not affect the results much, ranging from almost no effect at all for A. aurita to a slight drop in neutrality for the seawater samples, from R² = 0.7 at 1,000 reads/sample to R² = 0.6 at 50,000 reads/sample. Only for the communities associated with the sponge I. oros and two of the C. foliascens populations did read depth show a more pronounced effect; in both cases, neutrality dropped from R²≈ 0.6 at 1,000 reads/sample to R²≈ 0.45 when read depth was exceeding 10,000 reads/sample. Interestingly, an opposite trend was observed for the sediment samples, in which neutrality increased from R²≈ 0.5 at 1,000 reads/sample to R²≈ 0.7 at 10,000 reads/sample. The minimal to moderate changes in neutrality with increasing read depth for some datasets potentially reflect the influence of rare, non-neutral taxa, which are only detected with higher read depths. Shaded areas indicate 95% bootstrap confidence intervals. The data used to generate this figure are available in S1 Data.

https://doi.org/10.1371/journal.pbio.3000298.s009

(PDF)

S7 Fig. Overlap of taxa found in hosts and the environment.

Left: for two sponge species, there is only a very small overlap between the sponge microbiota and the taxa found in seawater. Right: for C. elegans, only a subset of the environmentally available microbes is found in the worms. The OTU abundance data used to generate this figure are available in S2–S10 Data. OTU, operational taxonomic unit.

https://doi.org/10.1371/journal.pbio.3000298.s010

(PDF)