Introduction

Hubbell’s1 neutral theory presented a null model for testing the mechanisms of species coexistence and biodiversity maintenance in ecological communities. It also offered a simple mechanistic explanation of the species abundance distributions (SAD), which has been extensively studied in community ecology since the 1940s, but still lacks a theoretical consensus for their interpretations2,3,4,5,6,7,8. Neutral theory contrasts with the traditional niche theory9,10,11,12,13,14 by assuming that the properties of individuals in a community are independent of their species identity. The traditional niche theory would assume that individuals of different species occupy different niches due to their species-specific properties and the niches are of limited similarities. According to neutral theory, local communities are built by random draws from regions pools driven by stochastic colonization, and the deterministic competitive interactions are insignificant in shaping the local community compositions because species are supposed to be competitively equivalent1,15,16. The neutral theory also assumes that abundant regional species possess a higher colonization probability and consequently the relative abundance of each species in a local community should mirror its abundance in the regional community at any given time1,16. Consequently, neutral processes driven community should be less influenced by the differences in environmental factors, and demonstrate less of a direct link between the environment and assembly processes17. The significance of Stephen Hubbell’s neutral theory of biodiversity2 is multifaceted, whether or not the null model fits to a specific community. For example, as demonstrated in this article, the traditional view that “everything is everywhere, but the environment selects” originally proposed by Baas-Becking18 is a well-known doctrine but receiving occasional challenges especially in recent years, indeed make great sense from the neutral theory perspective.

During the past decade, there has been extensive testing of the neutral theory, mostly with macro community ecology data19,20,21,22,23. McGill21 reviewed empirical tests of neutral theory; many of the tests used the neutral theory model as a null model and used another model (often lognormal distribution model) as the alternative model. However, in many cases, the alternative model was loosely defined, and hence the rejection of neutral theory could not automatically prove that a particular alternative niche-based model holds. In fact, the approach of using simplified assumptions, which are not required to be strictly accurate, about complex systems to improve understanding and to see what can be explained with the simplified models is a well-accepted practice in many fields of sciences24. Therefore, the primary value of the neutral theory is not whether or not it can adequately describe a particular dataset because the failure itself can inform us that processes other than the neutral processes prevail in shaping the structure of the community under investigation.

The neutral theory was developed and tested primarily in the context of ecological communities of plants and animals. In contrast to the many macro-ecological studies, there are relatively few applications of neutral theory in microbial ecology17,18,25,26,27,28,29,30,31,32,33,34,35,36,37,38. The scarcity is primarily due to the reality that majority of microbes are not cultivable, and their detections were extremely difficult until the invention of the metagenomic technology nearly a decade ago, which revolutionized the techniques for studying microbial community ecology39,40,41,42,43. When majority of the microbes could not be counted, or even detected, it would be very difficult if not impossible to collect sufficient data for testing ecological theory, including the neutral theory. Today the wide availability of metagenomic technology and the on-going Human Microbiome Project (HMP) offer us an unprecedented opportunity to apply and test major ecological theories for the studies of microbial communities43.

Limited tests of the neutral theory with microbial communities have produced mixed results in the existing literature. Earlier studies in environmental microbiology were mostly positive, supporting the neutral theory (e.g. refs 26, 27, 28), but more recent studies, especially those studies on the human and animal microbiomes seem to be mixed with more negative cases17,30,32,36,44 than positive cases35,37. Sloan et al.26 and Woodcock37 were among the first who explored the roles of immigration and chance in shaping prokaryote community structure. Woodcock28 performed the first comprehensive testing of the neutral theory in microbial ecology, and presented strong evidence that neutral community models fitted to the microbial communities in tree holes filled by rain water28,45 conjectured that the scale of regional communities in microbes might be different from that of macro-communities since the dispersal distance of microbes can be more extensive. In practice it is most likely that both niche and neutral mechanisms are in effect in microbial communities46,47,48,49,50, but niche effect is often more prevalent than neutral effect46,47. Of course, there were also the opposite cases, i.e., neutral effect is more significant (e.g. ref. 48).

The potential importance of neutral process in shaping the human microbiota has been conjectured in some perspective reviews43,51, but a handful of existing tests of the neutral theory with the human microbiome (including animal microbiome) only offer conflicting evidence. Jeraldo47 tested the role of neutral process in structuring the gut microbiome of three domesticated vertebrates, and they concluded that, although the species abundance patterns were seemingly well fit by the neutral theory, the theory couldn’t explain the evolutionary patterns in the genomic data (i.e., phylogenetic diversity). Furthermore, their analyses strongly supported the non-neutral (niche) role in shaping the animal microbiomes. Avershina et al.31, through a study of the faecal microbiota with 16S rRNA sequencing technology from a healthy cohort of 86 mothers and their children, observed a clear age-related colonization pattern shift in children and suggested that neutral processes are involved in shaping the gut microbiota. In contrast, there was not a similar shift in microbial composition during mother’s pregnancy. This finding is not consistent with previously mentioned finding from Jeraldo47 study on the gut microbiome of three domesticated vertebrates. Avershina et al.31 suggested that the inconsistence could be reconciled by the fact that niche selection should limit the phylotypes allowed in a given environment (the gut), whereas among the allowed phylotypes, neutral process could be in effect. Levy & Borenstein32 demonstrated the non-neutral assembly processes and complex co-occurrence patterns in the human gut microbiome. Guan & Ma44 tested the neutrality of human milk microbiome, and the test results were negative. Schmidt et al.17 found that deterministic processes appear to govern the assembly of fish microbiome, and stochastic colonization does not occur in their system. In the above four case studies17,32,44,47, the human or animal microbiome datasets failed to support the neutral theory.

Besides above mentioned gut microbiome study by Avershina et al.31, the following two case studies also presented supporting evidence to the neutral theory. Morris et al.35 applied the neutral theory to test whether dispersal from the mouth (source community) can adequately describe the observed microbial distribution in the lung. They also tried to identify microbial groups that significantly deviate from neutral community pattern, and they conjectured that these microbes might be advantageous in the lung habitat. They discovered that lung bacterial populations were similar to those in the oropharynx, but the lung does have some distinctive bacterial populations, not from contamination or dispersal from the mouth. Nevertheless, the functions of these lung habitat specific bacteria as well as the nature of the immune response to them are open for further research. Venkataraman et al.37 further utilized the neutral theory model to distinguish the microbes dispersing to the lungs from other body sites and atmosphere from the microbes indigenous to the lungs, assuming that the former (“neutrally distributed” species, i.e., consistent with dispersal and ecological drift) can be predicted with a neutral model and consequently the latter can be predicted from the difference between the actual presence and the former. They also suggest that the non-native microbial species should have a competitive advantage to the indigenous lung microbes. Their study reveals that the bacterial community composition of the healthy lung can be satisfactorily predicted with the neutral model, reflecting the overriding role of dispersal of microbes from the oral cavity in shaping the microbial community in healthy lungs. In contrast, it is not the case with the microbiome of diseased lungs where active selection is underway. Therefore, the bacterial neutral distribution is a distinguishing characteristic of the microbiome in healthy lungs. But, bacterial populations in diseased lungs appear to be under active selection. The ability to measure the relative importance of selection and neutral processes including dispersal in shaping and maintaining the healthy lung microbial community is a critical advance for understanding its influence on host health37. Ding & Schloss52 demonstrated that it is possible to distinguish different community types despite considerable intra-subject and inter-subject variation in the human microbiome, which has important implications for assessing disease risks associated with certain community types. If neutral theory can be utilized to distinguish healthy and diseased microbial communities, such as demonstrated by Venkataraman et al.37, its practical biomedical significance is self-evident.

More recently, advances have also been made in the theoretical exploration of the neutral process in microbial communities. Holmes et al.30 and Harris et al.36 reformulated Hubbell2 neutral theory model as a hierarchical Dirichlet multinomial mixture process to simulate the human gut microbiome. They concluded that functional niches mainly determine the human gut microbiome, and neutral assembly may only play a partial role in shaping the fine-scale gut microbial diversity—operating within the species occupying a specific function niche (i.e., those species playing a same metabolic role). They termed those cases “borderline” neutral patterns. They also discovered a negative correlation between body mass index (BMI) and immigration rates within the family Ruminococcaceae36, which offers a freshly new interpretation of the relationship between obesity and gut microbiome. Harris et al.36 not only revealed these important insights on gut microbiome but also developed a general computational procedure to efficiently fit multisite neutral theory models. We will further discuss this significant methodological advance in the discussion section.

O’Dwyer et al.33,34 tested the neutral theory from phylogenetic diversity perspective, rather than traditional species diversity perspective, with human microbiome data. Interestingly, they found a clear impact of metacommunity size (scale) on the phylogenetic diversity of body habitats relative to the null hypothesis. They also found that whole microbiome diversity for a given subject is typically much lower than a random sample from the metacommunity, which is complicated by a wide range of different behaviors for the distinct habitats within that subject. Overall, they concluded that there are significant variations in diversity across habitats relative to the prediction from the neutral theory. Zeng et al.38 developed an agent-based architecture that encapsulated the neutral model of community diversity with added dimensions of an evolutionary timescale and a genealogy of hosts. This is essentially a simulation system that can be utilized to test various assumptions on community assembly including neutral theory assumptions and its advantages include the consideration of the parental contribution effects on community assembly process as well as how the process operates over evolutionary time.

The objective of this study is two-fold: First, we conduct a comprehensive testing, with extensive HMP datasets covering all five major body sites (gut, oral, skin, vaginal, and nasal) of the human microbiome, including 7437 communities from 18 locations of 242 individuals, and with rigorous statistical simulation tests to answer the question—is the community assembly in the human microbiome neutral? Second, we discuss the implications of the neutral theory to understanding the human microbiome diversity and explore the underlying mechanisms if the data fail to support the neutral theory.

Materials and Methods

The 16s rRNA sequence datasets of human microbiome

The datasets we use to test the neutral theory were obtained from HMP data center (www.hmpdacc.org). Specifically, we selected the 16s rRNA sequence datasets of 18 sites of 5 locations sampled from 242 subjects42 (termed 240-healthy-subjects HMP datasets, hereafter). Those body sites and locations are: airways (anterior nares), gut (stool), oral (attached keratinized gingival, buccal mucosa, hard palate, Palatine tonsils, saliva, subgingival plaque, supragingival plaque, throat, tongue dorsum), skin (left antecubital fossa, left retroauricular crease, right antecubital fossa, right retroauricular crease), and urogenital (mid vagina, posterior fornix, vaginal introitus). We used the datasets of both V1-V3 region and V3-V5 region in our analyses, and totally 7437 community samples (V1-V3, 2855 samples; V3-V5, 4582 samples, respectively) were sequenced, and two collections of OTU tables (corresponding to V1-V3, and V3-V5 regions, respectively) were obtained from the 16s-rRNA mothur software pipeline. Each sample corresponds to one row in the OTU tables, and is used to fit the neutral theory model. The OTU tables were computed with a cutoff 3% dissimilarity based on RDP database and they are available for download at www.hmpdacc.org. Those OTU tables contain species abundances (OTU reads) within each community sample.

The computational procedures

The essential aspects of the neutral theory can be summarized as: (i) it assumes that interacting species are equivalent on an individual ‘per capita’ basis, (ii) it is an individual-based stochastic dynamic theory, (iii) it is a sampling theory, and finally (iv) it is a dispersal-assembly theory53. The theory reveals the significance of dispersal limitation, speciation and ecological drift in community assembly and maintenance. As a null model, it offers us an apparatus to evaluate the role of adaptation and natural selection in the context of evolutionary community ecology54. The classic neutral theory model describes a local community containing J individuals. One of these individuals, chosen at random, dies and is replaced at every time step. The replacement can be either an offspring of another randomly chosen individual from the local community, occurring with probability 1-m, or an offspring of a randomly chosen individual from the metacommunity, with probability m. The parameter m is considered as a measure of dispersal limitation. A problem with m is that it does not translate into dispersal limitation in the most logical way; for example, a small local community may involve a much smaller flux of immigrants for the same value of m. For this reason, an alternative parameter I is often used instead of m55.

We use two sampling formulae for testing the neutral theory: one was proposed by Ewens56, and another by Etienne57 that was also inspired by Ewens formula. Both are accurate likelihoods for different scenarios, and Etienne’s one is for the best-known model of Hubbell’s2 neutral theory. Etienne’ s formula adds dispersal limitation to Ewens formula, and we can infer that dispersal limitation plays a significant role in neutral community assembly if Etienne’s formula performs better than Ewens formula. In literature, the comparisons between Ewens and Etienne and sampling formulae showed that the latter outperformed the former when the dispersal limitation plays a significant role in community assembly57.

We will first compare the log-likelihoods computed with Ewens formula and Etienne formula to determine which of them is better suited for the human microbiome datasets. We then compare the log-likelihoods of the observed community (sample) and neutral-theoretic community to test whether or not the neutral theory model fits to the observed community, based on either Ewens56 or Etienne57, whichever performed better in the previously described first step. In addition, we also hope that the possible difference between Ewens56 or Etienne57 and log-likelihoods will shed light on the effect of dispersal limitation.

Ewens sampling formula was proposed to describe the probability distribution of alleles in genes in the context of molecular neutral evolution56,58,59,60. Hubbell2 applied it to compute the likelihood of a given ecological community to satisfy the prediction of the neutral theory:

where, J is the total number of individuals in the community, S is the total number of species, θ is the biodiversity parameter of the sampling formula, niis the abundance of species i, φa is the number of species with abundance a. Hubbell2 had realized the important role of dispersal limitation but the sampling formula he originally used, i.e., Ewens formula, did not consider the effect of dispersal limitation (i.e., immigration rate m = 1).

Etienne57 proposed a new sampling formula with dispersal limitation (m < 1),

The resulting Etienne57 sampling formula is

where K(D, A) is defined as

Etienne57 compared the two sampling formulae using Barro Colorado (BCI) tree dataset61 and Caño Maraca (CM) fish dataset62, and the results showed that Etienne sampling formula led to a significantly better fit to the datasets than Ewens sampling formula did. Here we compare the two formulae with HMP datasets to detect the potential effect of dispersal limitation. The following log-likelihood ratio test is used to compare the results from both the formulae:

where, L0 is the null model and L1 is the alternative model, D is the deviation that is twice the difference between the log-likelihoods of the two formulae. The p-value computed follows a X2-distribution with the degree of freedom being one.

The method we use to test the community neutrality is Etienne’s63Exact neutrality test.” The exact neutrality test is based on the sequential construction approach, and an obvious advantage of this method is that no alternative model is needed in the approach. The test is conducted as follows: Firstly, maximum likelihood estimation (MLE) is used to estimate the model parameters from the observed data. Secondly, the set of model parameters [θ, I] estimated with MLE and parameter (J, the total number of individuals in the community) can be used to generate any number of artificial datasets (D). In our case of HMP data, 100 artificial datasets for each community sample are generated. Thirdly, we use Etienne’s formula to calculate the likelihood of artificial datasets, and compare the average likelihood of artificial data with that of the observed data.

The sampling formula [(1) or (3)] gives the probability (PS) of any dataset in D (|D| = 100 in our case) satisfying the neutral model. The sampling formula also gives the probability P0 of real observed dataset satisfying the neutral model. If the probability computed from the real dataset is significantly smaller than the average computed from the artificial datasets, it is unlikely that the observed community pattern was structured by the neutral process. If the probability from observed community is close to the values computed from the artificial datasets, then the observed species abundance distribution (SAD) is consistent with neutrality, and we conclude that the neutrality of the community under testing cannot be rejected.

Results and Discussion

In the following, we report our results in two parts, corresponding to (i) the comparison of Ewens and Etienne sampling formulae, which determines which formula will be used to test the community neutrality in the next step; (ii) the comparison of the log-likelihoods of theoretical and observed communities to determine whether or not an observed community is neutral.

Comparison of Ewens and Etienne sampling formulae

For each human microbial community sample, we compute the fundamental biodiversity parameter (θ), the immigration probability (m), and the likelihood of seeing the dataset given these parameter values, with Ewens and Etienne sampling formulae respectively. For each community sample, the 16s rRNA sequence data of both V1-V3 and V3-V5 variable regions were separately utilized to test the neutral theory respectively. Tables 1 and 2 are the summary results of the comparisons between Ewens and Etienne sampling formulae for region V1-V3 and V3-V5, respectively. Detailed comparative results of the two formulae are provided as supplementary Tables S1 and S2, in which parameters and statistics such as the biodiversity parameter (θ), the immigration probability (m), log-likelihood, and likelihood ratio are provided to demonstrate the suitability of the both formulae.

Table 1 The comparison of Ewens and Etienne sampling formulae with V1-V3 region data.
Table 2 The comparison of Ewens and Etienne sampling formulae with V3-V5 region data.

In Tables 1 and 1, the first column lists the five sampling location, and the second column lists the 18 sites that are distributed over the five sampling locations. Column 3 is the total number of communities for each site. Column 4 is the number of communities for which both Ewens and Etienne formulae made no significant difference (NSD). The results from Tables 1 and 1 show that in both V1-V3 and V3-V5 regions, there are little differences between Ewens and Etienne sampling formulae in terms of the log-likelihood (p > 0.05). The percentages of NSD communities in both Tables 1 and 1 are close to 100%, i.e., both formulae made no differences in almost all human microbial communities we tested.

Etienne’s neutrality exact test

One hundred (100) artificial communities were generated to simulate the theoretical neutral community corresponding to each observed human microbial community. As mentioned previously, we have 2855 communities represented by V1-V3 data, and 4582 communities represented by V3-V5 data. The parameters (m, θ, J, S) of the neutral theory models for those human microbial community that have passed the likelihood ratio test are listed in Tables 3 and 3, for V1-V3 and V3-V5, respectively. It is noted that Etienne formula63 was utilized for computing the results in both Tables 3 and 3. In fact, Ewens formula could be used either since both the formulae made no difference.

Table 3 The 26 human microbial communities that passed the neutrality exact tests with Etienne sampling formula (V1-V3 region)*.
Table 4 The 23 human microbial communities that passed the neutrality exact tests with Etienne sampling formula (V3-V5 region)*.

Table 3 lists the parameters of 26 human microbial communities, based on V1-V3 region data, which have passed the neutrality exact test, and the results of the other 2828 communities that failed to pass the neutrality exact tests, were omitted in Table 3 (but provided in the Supplementary Table S3). It also shows that 99.1% human microbial communities based on the 16s rRNA sequence data of V3-V5 region do not satisfy the neutral community model.

Similarly, Table 4 contains the parameters of 23 human microbial communities, based on V3-V5 region data, which have passed the neutrality exact test, and the results of the other 4557 communities that failed to pass the neutrality exact tests, were omitted in Table 4 (but provided in the Supplementary Table S4). It also shows that 99.5% human microbial communities based on the 16s rRNA sequence data of V3-V5 region do not satisfy the neutral community model.

Note that among the 26 communities that passed the neutrality exact tests based on V1-V3, 21 are urological microbiome. In other words, the majority (80%) of the neutral human microbial communities are urological microbiome. Among the 23 communities that passed the neutrality exact test based on V3-V5, 14 are skin microbiome. In other words, the majority (60%) of the neutral human microbial communities are skin microbiome.

In summary, our results show that less than 1% (49 out of 7437 community samples) satisfied the neutral theory model according to the tests with Ewens56 and Etienne57 models. We therefore conclude that the human microbiome is not neutral in general. Figure 1 shows the graphs of four samples that successfully fit to the neutral theory model of Etienne’s63 neutrality exact test.

Figure 1: The rank abundance curves of the four samples, out of the 49 samples that were successfully fitted to Etienne’s neutral model.
figure 1

The four samples represent the four different sites from four different individual subjects: throat (Subject ID: 700110166), left antecubital fossa (Subject ID: 700102129), right antecubital fossa (Subject ID: 700106174), and stool (Subject ID: 700106979). In the figure, the solid red lines are observed data and the black dash lines are artificial (simulated) datasets. The X-axis is the species rank order in abundance and Y-axis is the abundance of each species in natural logarithm.

Discussion

Bell64 proposed a method termed distance-decay dispersal for testing the dispersal limitation, which analyzes the relationship between the distance and community dissimilarity. If there is dispersal limitation effect within the community, then community dissimilarity should increase with the increase of distance, and the relationship between distance and dissimilarity should be linear64,65. Otherwise, the relationship should be non-linear. By simulation studies with two simple theoretical community models, Lotka-Volterra stochastic dynamics model and presence-absence model, Fisher & Mehta66 discovered that there is a phase transition in diverse ecological communities between a selection-driven regime (the niche phase) and a drift-dominated regime (the neutral phase), similar to the phase transitions occurring in H2O (i.e., solid ice, liquid water, and gas steam). Their simulation study suggests that the niche phase is more likely in communities with large populations sizes and relative constant environments, and the neutral phase with small population sizes and fluctuating environments. The authors admitted that some simplifications, which may be unrealistic for natural ecological communities, were introduced in their study: one example is that the community is well mixed with purely competitive interactions. In reality, most natural communities are patchy, hierarchical, of complex spatial structures, and the effects of dispersal mutation, and mutualism can hardly be ignored, especially in microbial communities. Indeed, if the prediction from their simulation is taken literally, the tropical forests, which should have relatively stable environment, should be less likely to be neutral. However, it is well known that Hubbell2 neutral theory was strongly inspired by his studies of tropical forests. Nevertheless, we concur with Fisher & Mehta66 that the presence of a niche–neutral phase transition is likely robust to the additional model modifications with more realistic assumptions. We also fully agree with them that more realistic modifications are needed to develop more quantitative models for natural communities.

Inspired by Fisher & Mehta66 study, we propose the following hypothesis to explain the failure of neutral theory model in describing the human microbial communities. The famous microbial biogeographical doctrine first proposed by Baas-Becking18, “Everything is everywhere, but the environment selects” offers us another piece of key supporting evidence, besides Fisher & Mehta’s66 finding. We conjecture that the balance between dispersal and selection may tip the shift (phase transition) between the neutral and niche regimes. While dispersal favors neutral processes, selection obviously weighs against the effect of neutral processes and in favor of the niche phase.

Phylogenetic diversity has increasingly become an important tool to quantify community diversity thanks to the rapid accumulations of the metagenomic sequencing data67,68. O’Dwyer et al.33,34 introduced a framework that uses the phylogenetic diversity for testing the neutral theory. They compared patterns of phylogenetic diversity across multiple bacterial community samples from different habitats with the evolutionary trees generated using theoretical models of biodiversity. They obtained two major findings: the first finding is that on coarse scales the backbone of the empirical trees is simple, robust and consistent across habitats, although they are idiosyncratic on finer scales. While their first finding still supports the primary principle of the neutral theory—that selective differences are irrelevant for predicting large-scale patterns, and therefore the neutral theory is likely to predict the biodiversity patterns over coarse scales, their second finding reject the neutral theory model per se. They found that the existing neutral models could not explain observed patterns in microbial phylogenetic diversity, and instead, another family of Λ-coalescents models offers a qualitatively better description of both the scaling and topology of empirical trees. Their overall conclusion is that while the principle of the neutral theory is useful as a backbone for characterizing the microbial phylogenetic diversity, new generation of neutral models such as the family of Λ-coalescents models are needed to implement the utilization of the backbone for revealing the ecological and evolutionary mechanisms of microbial diversity. While we fully agree with O’Dwyer et al.33,34 that testing of the neutral theory with phylogenetic diversity makes great sense, especially with new type of neutral models, we suggest that different metrics for diversity may also make difference. Two biodiversity metrics should be particularly worthy of pursuing in future testing of the neutral theory. One is the UniFrac by Lozupone & Knight67, and another is the Hill numbers. The latter has increasingly been recognized as the most appropriate alpha-diversity metrics, and its multiplicative partition for beta-diversity is found to have excellent statistical and ecological properties advantageous over other diversity indexes68,69,70.

Metagenomics technology opens unprecedented opportunities to test ecological theories such as the neutral theory of biodiversity by producing gigantic amount of molecular sequencing data often from many community samples. This nevertheless also presents significant computational challenges. Some of the challenges can be dealt with standard bioinformatics software tools such as QIIME and Mothur71, but a far more significant challenging is to maximally take advantage of the big data of metagenomic sequences. One such challenge in the case of neutral theory modeling is the extreme computational complexity involved in multisite neutrality testing, which is computationally intractable when the number of sites (local communities) are more than a few36. A recent major methodological advance by Harris et al.36 offers a timely approach to the problem. The approach approximates a large class of neutral models with the hierarchical Dirichlet process by developing a highly efficient Bayesian fitting strategy for the multisite neutrality-testing problem. Besides being able to handle large datasets in a reasonable amount of time for multi-site neutrality test, an additional benefit is to generate full posterior distribution over the parameters, rather than obtaining just a maximum likelihood prediction. The approach also reconstructs the metacommunity distribution, which makes it possible to divide the problem of neutrality test into two parts. First, from the full neutral models with fitted parameters, one can generate samples and compare their likelihood with that of the observed samples to test for neutrality. Second, one can generate samples from the observed metacommunity and test for the neutrality of local community alone. Therefore, it is possible to test for neutrality at both local and metacommunity level with Harris et al.36 hierarchical Bayesian modeling approach. From a broader perspective, it is suggested that to use hierarchical Dirichlet process as an ecological null model, possibly extended to niche-neutral hybrid model and playing a more significant role in the community ecology36.

Additional Information

How to cite this article: Li, L. and Ma, Z. (S.) Testing the Neutral Theory of Biodiversity with Human Microbiome Datasets. Sci. Rep. 6, 31448; doi: 10.1038/srep31448 (2016).