Testing the Neutral Theory of Biodiversity with Human Microbiome Datasets

The human microbiome project (HMP) has made it possible to test important ecological theories for arguably the most important ecosystem to human health—the human microbiome. Existing limited number of studies have reported conflicting evidence in the case of the neutral theory; the present study aims to comprehensively test the neutral theory with extensive HMP datasets covering all five major body sites inhabited by the human microbiome. Utilizing 7437 datasets of bacterial community samples, we discovered that only 49 communities (less than 1%) satisfied the neutral theory, and concluded that human microbial communities are not neutral in general. The 49 positive cases, although only a tiny minority, do demonstrate the existence of neutral processes. We realize that the traditional doctrine of microbial biogeography “Everything is everywhere, but the environment selects” first proposed by Baas-Becking resolves the apparent contradiction. The first part of Baas-Becking doctrine states that microbes are not dispersal-limited and therefore are neutral prone, and the second part reiterates that the freely dispersed microbes must endure selection by the environment. Therefore, in most cases, it is the host environment that ultimately shapes the community assembly and tip the human microbiome to niche regime.

of sciences 24 . 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 ecology 17,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 ecology [39][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 communities 43 .
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 cases 17,30,32,36,44 than positive cases 35,37 . Sloan et al. 26 and Woodcock 38 were among the first who explored the roles of immigration and chance in shaping prokaryote community structure. Woodcock 28 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 water 28,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 communities [46][47][48][49][50] , but niche effect is often more prevalent than neutral effect 46,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 reviews 43,51 , but a handful of existing tests of the neutral theory with the human microbiome (including animal microbiome) only offer conflicting evidence. Jeraldo 47 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 Jeraldo 47 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 & Borenstein 32 demonstrated the non-neutral assembly processes and complex co-occurrence patterns in the human gut microbiome. Guan & Ma 44 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 studies 17,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 health 37 . Ding & Schloss 50 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 Hubbell 2 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 Ruminococcaceae 36 , 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 subjects 42 (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 theory 53 . 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 ecology 54 . 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 m 55 .
We use two sampling formulae for testing the neutral theory: one was proposed by Ewens 56 , and another by Etienne 57 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's 2 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 assembly 57 .
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 Ewens 56 or Etienne 57 , whichever performed better in the previously described first step. In addition, we also hope that the possible difference between Ewens 56 or Etienne 57 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 evolution 56,58-60 . Hubbell 2 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, n i is the abundance of species i, ϕ a is the number of species with abundance a. Hubbell 2 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). Etienne 57 proposed a new sampling formula with dispersal limitation (m < 1), The resulting Etienne 57 sampling formula is Etienne 57 compared the two sampling formulae using Barro Colorado (BCI) tree dataset 61 and Caño Maraca (CM) fish dataset 62 , 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, L 0 is the null model and L 1 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 X 2 -distribution with the degree of freedom being one. The method we use to test the community neutrality is Etienne's 63 "Exact 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 (P S ) of any dataset in D (|D| = 100 in our case) satisfying the neutral model. The sampling formula also gives the probability P 0 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 Scientific RepoRts | 6:31448 | DOI: 10.1038/srep31448 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.
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  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 formula 63 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 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 Ewens 56 and Etienne 57 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's 63 neutrality exact test.   Discussion Bell 64 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 linear 64,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 & Mehta 66 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 H 2 O (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 Hubbell 2 neutral theory was strongly inspired by his studies of tropical forests. Nevertheless, we concur with Fisher & Mehta 66 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 & Mehta 66 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-Becking 18 , "Everything is everywhere, but the environment selects" offers us another piece of key supporting evidence, besides Fisher & Mehta's 66 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 data 67,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 & Knight 67 , 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 indexes [68][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 Mothur 71 , 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 few 36 . 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 ecology 36 .