The relationship between human vaginal microbiome (HVM) and women’s health has been investigated since the 1980s, when clinical microbiologists had postulated that the diversity and possibly stability of vaginal microbiome are involved in the occurrence/recurrence of bacterial vaginosis (BV) (e.g., Sobel, 1999; Fredricks et al., 2005; Fredricks, 2011; Ma et al., 2012). Those studies are among the earliest ecological approaches to diseases now often referred to as the human MADs (microbiome-associated diseases) with an ever more rapidly growing list including BV, IBD (inflammatory bowel disease), periodontitis, cystic fibrosis (CF), psoriasis and many others (Lynch and Pedersen, 2016; Knight et al., 2017; Young, 2017; Gilbert et al., 2018). The metagenomics technique and the launch of the human microbiome project (HMP) and MetaHIT (metagenome of human intestinal tract) have revolutionized the investigation of the human microbiome and associated diseases during the last decade or so. Nevertheless, many questions in the field are still open and new more complex questions are being raised. In the case of BV and vaginal microbiome, as described by Fredricks (2011) who borrowed Winston Churchill’s words for a very different topic, “BV remains a riddle, wrapped in a mystery, and inside an enigma.” A recent characterization “that BV is not a single entity, but a syndrome linked to various community types that cause somewhat similar physiological symptoms.” by Ma et al. (2012) reflects the state-of-the-art understanding of BV etiology. Obviously, although the importance of vaginal microbiome ecology in BV etiology is repeatedly confirmed, the mechanistic relationship between BV and HVM is far from clear. A pair of questions of fundamental importance: what are the underlying mechanisms driving the dynamics of HVM and what are their implications to the occurrence/recurrence of BV, are still largely unanswered.

Addressing the question of community assembly and diversity maintenance, the essential ingredients of community structure and dynamics, has attracted extensive attention and also led to vigorous debate (Alonso et al., 2006; McGill et al., 2006; Chisholm and Pacala, 2010; Rosindell et al., 2011). Two leading and competing theories in this field have been the traditional niche theory with a history back to the 1910s (Grinnell, 1917; Hutchinson, 1957; Holt, 2009) and more recent neutral theory (Hubbell, 2001). Both theories were invented to explain a familiar phenomenon on the earth, which was described by Darwin (1859) in the last paragraph of his “On the Origin of Species” as “It is interesting to contemplate a tangled bank, clothed with many plants of many kinds, with birds singing on the bushes, with various insects flitting about, and with worms crawling through the damp earth, and to reflect that these elaborately constructed forms, so different from each other, and dependent upon each other in so complex a manner, have all been produced by laws acting around us.” In modern ecological terminology, entangled bank is essentially the concept of ecological community. Darwin was wondering how diverse lives (species) could coexist and form a beautiful entangled bank, while his theory stipulated the universal struggle for life as a consequence to natural selection. The classic niche theory assumed that each species has its own niche in which its individuals are adapted to live and prosper, and the entangled bank consists of many different niches suitable for many different species. In terms of niche theory, deterministic traits a species possess or selective niche forces play critical roles in driving the assembly of an ecological community as well as the maintenance of diversity after the community is established.

In the late 1990s, Hubbell (2001) challenged the traditional niche view by proposing the unified neutral theory of biodiversity and biogeography (UNTB). Different from traditional niche theory, the UNTB was formulated as a probability distribution model, which can be fitted with the species abundance distribution data (the number of each species in a community), obtainable by sampling ecological communities, and rigorously tested statistically. The theory assumes that the individuals of all species in a community are demographically equivalent, but their birth/death rates are stochastic, which means birth-death, migration, and speciation are all random events. Consequently, random drift and dispersal play critical roles in driving community assembly and diversity maintenance. Some researchers argued that the concept of species equivalence is “flawed” given the existence of niche differences and competitive asymmetries among species. Nevertheless, the stochasticity in species demography (particularly of single-cell microbes) is also a biological reality and its role may not be ignored in many communities. In reality, both deterministic niche forces and stochastic neutral forces may be in effect in setting the rules of community assembly and diversity maintenance, and it may be the hybrid effects that shape the community dynamics. For this reason, in the last decade and so, several hybrid models that integrate neutral and niche effects have been developed (e.g., Tilman, 2004; Ofiteru et al., 2010; Stokes and Archer, 2010; Jeraldo et al., 2012; Pigolotti and Cencini, 2013; Tang and Zhou, 2013; Fisher and Mehta, 2014; Kalyuzhny et al., 2014a, b, 2015; Matthews and Whittaker, 2014; Noble and Fagan, 2015). As to the debates on the usefulness and validity of the UNTB, using an analogy, in modern statistics (especially in biostatistics), it has been widely recognized that many datasets do not follow the Gaussian distribution (the normal distribution); nevertheless, few statisticians would question the foundational role of the Gaussian distribution, not to mention its validity. Similarly, the merits and unique advantage of UNTB as a null model for testing the significance of stochastic drift and dispersal have been firmly established and widely applied in the community ecology of plants and animals.

In the present study, we use a pair of models, the first a multi-site neutral model (Harris et al., 2017) and the second, a niche-neutral hybrid model (Tang and Zhou, 2013), to evaluate the relative significance of neutral and niche effects in shaping the dynamics of HVM. We further investigate the difference in the neutral-niche continuum between BV patients and healthy women. Our approach is different from most existing applications of neutral or hybrid models in the following three aspects.

First, most existing neutral or niche-neutral hybrid models use spatially implicit/explicit community/metacommunity data, whereas we use longitudinal (time-series) sampling of the community/metacommunity. In spatially explicit models, the metacommunity consists of multiple local communities, which are connected with each other through dispersal and migration. In temporal (time-series) models, the metacommunity consists of a series of “snapshots” of the same community at different time points, i.e., the time-series data obtained from sampling the vaginal microbiome of a subject at different time points in this study. Indeed, previously, Kalyuzhny et al. (2014a, b, 2015) used time-series data to perform dynamic analysis of the niches versus neutrality and they termed the analysis as a generalized neutral theory for explaining the static and dynamic properties of ecological communities. A reason we did not adopt their models is that the models we use in this study, as explained below, are truly multi-site mechanistically, which are mapped to the time-series points in our study.

Second, we use a truly multi-site neutral (MSN) model of UNTB, which was developed by Harris et al. (2017) to overcome the severe computational limitation of existing neutral theory models when the number of local communities is large and the migration rates among the local communities are different (Etienne, 2007, 2009a,b). The core technique Harris et al. (2017) developed was to approximate the multi-site UNTB model with the hierarchical Dirichlet process (HDP) and use an efficient Bayesian machine-learning algorithm. With their approach, fitting even the largest dataset can be performed in a reasonable amount of time. This important computational advance enables us to build a UNTB model for each subject by utilizing the time series sampling of her vaginal microbiome. This capability is of significant practical importance given the established connection between BV and the diversity of the vaginal microbiome, in particular, a long-standing puzzle in BV etiology – the rise of species diversity associated with BV (e.g., Sobel, 1999; Fredricks et al., 2005; Fredricks, 2011; Ma et al., 2012; Ma and Ellison, 2018, 2019).

Third, we also apply the niche-neutral hybrid (NNH) model by Tang and Zhou (2013) to further assess the neutral-niche hybrid effects in shaping the dynamics of HVM diversity. A major reason we prefer this hybrid model to other existing hybrid models (e.g., Tilman, 2004; Ofiteru et al., 2010; Stokes and Archer, 2010; Jeraldo et al., 2012; Pigolotti and Cencini, 2013; Tang and Zhou, 2013; Fisher and Mehta, 2014; Kalyuzhny et al., 2014a, b, 2015; Matthews and Whittaker, 2014; Noble and Fagan, 2015) is because both MSN and NNH use exactly the same data collection methods – either multi-site or multi-time-point sampling. The only essential difference between the neutral-niche hybrid model (NNH) and multi-site neutral model (MSN) is the assumption that niche differences exist among local communities in NNH, while the MSN assumes no niche differentiation. In a time-series setting, the NNH model can tell us whether deterministic forces (similar to habitat selection in a spatial setting) such as when changes in the host’s physiology significantly influence the dynamics of vaginal microbiome diversity over time.

In summary, by building and testing the MSN and NNH models for each subject, we are able to evaluate the relative importance of stochastic forces (e.g., neutral dispersal, drift, and stochastic diversification) vs. deterministic forces (e.g., microbial interactions, host genome and physiology, menses, etc.) in shaping the dynamics of community diversity. Furthermore, if we treat BV or health status as part of the host physiology, testing the MSN/NNH models can reveal the impact of BV on the dynamics of the HVM diversity (assuming that diversity change is the consequence of BV), or reveal the diversity changes that induce BV (assuming that diversity change is the cause of BV). Regardless of the causal assumption, our approach offers a useful tool for evaluating the mechanisms (niche vs. neutral) of the dynamics of HVM diversity as well as the factors affecting the balances between different mechanisms. We demonstrate our approach (see Figure 1) by using the datasets (see Table 1) from three separate longitudinal studies on the HVM, including 79 subjects sampled at 2,733 time points (Gajer et al., 2012; Ravel et al., 2013; Romero et al., 2014).

With traditional neutral theory of biodiversity, the spatially explicit or implicit model describing the metacommunity consisting of multiple local communities is the most frequently used metacommunity model for testing neutrality (Hubbell, 2001; Rosindell et al., 2011, 2012; Ma, 2020). The use of longitudinal community/meta-community samples to perform dynamic analysis of the niches vs. neutrality and to further assess and interpret the community static and dynamic properties by harnessing the neutral theory has been few but can be equally effective (Kalyuzhny et al., 2014a, b, 2015). The integrated modeling with MSN/NNH in previous sections demonstrated another approach to generalizing the neutral-theoretic analysis to temporal meta-communities. Furthermore, we take advantages of a recent advance in computational statistics made by Harris et al. (2017) HDP-MSN machine learning algorithm. The HDP-MSN overcomes a significant computational bottleneck that existed in estimating the migration rates (m) when the number of local communities is large, which prevented large-scale testing of the UNTB with truly multi-site datasets. Nevertheless, truly multi-site datasets are scarce, especially in the studies of the human microbiome, where community samples are usually taken from unrelated individuals, and therefore dispersal (migration) among individuals is unlikely to occur on ecological timescales. In this study, we use the time-series sampling data in place of spatial sampling data. That is, the vaginal microbial community of each subject was sampled in patients at varying numbers of time points (6–60, see Table 1). By using time-series data with the MSN and NNH models, one can effectively evaluate the levels of stochastic neutral forces and deterministic niche forces in driving the community dynamics. In the case of time-series data, stochastic neutral forces may include stochastic fluctuations in demography (in the birth-death processes of bacterial cell divisions and deaths), which is the analog to ecological drift in neutral theory. The deterministic forces in time series data can include diversity- or dominance-dependent regulatory forces for community stability (dynamics) (Ma and Ellison, 2018, 2019).

As explained in the previous section of results, the results from the integrated niche-neutral hybrid analyses under a default P-value = 0.05 with MSN and NNH models in this study seem to suggest that neutral drifts play a dominant role in driving the community dynamics of the HVM, while the deterministic niche differentiation is moderate (approximately 17% in terms of passing NNH test). As further elaborated below, the assessment of the relative significance of neutral vs. niche may be strongly influenced by the model-choice (MSN or NNH) and P-value thresholds. It should also be reiterated that the conclusions obtained from this study are from analyzing the temporal dynamics data of the HVM, rather than from analyzing spatial metacommunity samples as usually performed in community ecology. In other words, the local communities in our analyses are simply the snapshots of an individual woman’s vaginal microbiome dynamics. Therefore, niche differentiations are also “temporal differentiations,” which might be relatively weak due to the nature of longitudinal observations. Future studies with “orthodox” spatial metacommunity samples should shed a more comprehensive picture on the community assembly and diversity maintenance of the HVMs. In the remainder of this article, we discuss two uncertainties regarding the test of the neutral and/or niche-neutral hybrid models.

First, it is well known that a significant challenge in investigating the mechanisms of community assembly or distinguishing the neutral from niche effects is that multiple independents models with possibly different ecological assumptions about the mechanisms may produce similar goodness-of-fittings to the same datasets, which is termed the lack of full model specificity in previous sections. This can make the inferences of definite mechanisms from different models difficult since the mapping from assumptions to mechanisms may not be one-to-one. Table 6 listed the breakups of successful fittings of the MSN/NNH models with P-value thresholds of 0.05–0.95, classified as four groups including successfully fitted to “MSN-only,” “NNH-only,” “both MSN & NNH,” and “neither MSN nor NNH.” Here, we first discuss the breakups when P-value threshold for testing the MSN/NNH model is set to 0.05, and the other thresholds are discussed shortly below. As exhibited in Table 6, overall, the “MSN-only group” (fitted to MSN uniquely) takes about 68% (ranged between 47 and 100%) of all cases and “both MSN & NNH” group (non-unique fittings) takes about 29%, (ranged between 0 and 53%) and in less than 3% cases (2 out of 79 individuals) neither MSN nor NNH model was fitted successfully. Therefore, in the majority (68% or 54 out of 79 individuals), the MSN model was able to uniquely interpret the neutral dynamics of the HVMs, when the P-value threshold for testing MSN/NNH was set to 0.05.

Second, another dilemma that may lead to uncertainty in testing the neutral or niche-neutral hybrid models is the choice of P-value threshold. Traditionally, the P-value was set to 0.05 in testing the neutral theory; when P > 0.05, the null hypothesis or model (satisfying the MSN or NNH model) cannot be rejected. In other words, when P > 0.05, the observed community is considered indistinguishable from what the theoretical model predicts. In previous sections, P-value = 0.05 was termed traditional “default.” However, one may set P-value to other threshold values. The higher the P-value is, the more likely (the higher likelihood) that the community is consistent with the model prediction. That is, when the P-value is set to higher threshold values, it is more difficult to reject the null model. In terms of the neutrality test based on the MSN model, it implies that accepting neutral hypothesis is more reliable (conservative). In terms of the NNH model, it implies that accepting non-neutrality (niche differentiation) is more reliable (conservative). Consequently, when larger P-value thresholds are adopted, the confidence (reliability) to accepting the null model (MSN or NNH) is raised and the confidence to reject the null model (MSM or NNH) is lowered. Table 4 listed the passing percentages (strictly, should be stated as percentage indistinguishable from model prediction) from testing MSN/NNH when P-value was set to 0.05, 0.5, 0.9, and 0.95, respectively. Obviously, as shown in Table 4, higher P-values correspond to a lower passing percentage of MSN-neutrality tests. When the P-value threshold was raised to 0.95, the passing percentage of MSN-neutrality test declined to approximately 50%, while the percentage was 97.3% when the P-value = 0.05. Raising the P-value threshold from 0.05 to 0.95 is a rather dramatic increase of the confidence level for not rejecting the null neutral model (or accepting the null model), still nearly half the metacommunities (microbiomes of individuals) passed the MSN model, suggesting that the neutral drifts indeed play a significant role in shaping the dynamics of the HVM.

An interesting observation is that, when the P-value threshold was set to default 0.05 (Table 4), the passing percentages of MSM testing were not significantly different between different treatments (see Table 5 for Fisher’s exact test). However, when the P-value threshold was set to 0.95 (Figure 5), the differences between the BV group (including ABV and SBV) and healthy groups (HEA1, HEA2, and pregnancy) were significant (13.3–16.7% vs. 68.2–75%). The healthy groups exhibited significantly more neutral communities than BV groups (Table 5). This result is actually puzzling. A traditional view has been that the vaginal microbiomes associated with BV usually have higher community diversity than the healthy counterparts, possibly due to the loss of dominant species such as Lactobacillus (e.g., Ma et al., 2012). Since communities with dominant species tend to contain more asymmetric interactions, they could be less likely neutral. Therefore, this puzzling result appears to contradict with the traditional view. Of course, the relationship between dominant species and non-neutrality is not necessarily positive; we hope that future more mechanistic studies will resolve the apparent inconsistency.

The adoption of different P-value thresholds may also affect the previously discussed specificity (uniqueness) of MSN/NNH model fittings. Table 6 shows the breakups of four categories (MSN-only, NNH-only, both MSN & NNH, neither MSN nor NNH) under different P-values ranging from P-value = 0.05 to 0.95. It appears that the model specificity of MSN seems to increase with the increase of P-value threshold adopted. This result should be expected since the increased P-value should raise the confidence level for accepting the null (MSN or NNH) model, and the passing percentages judged with higher confidence would decline accordingly.

Finally, the two previously discussed uncertainties associated with testing the neutral or niche-neutral hybrid models can be interdependent and proper resolving them requires both ecological science and statistical art. There may not be a perfect solution for resolving those issues due to both the ecological complexity and the limitation of frequentist statistical approaches. The art lies in balancing the trade-off between reliability (confidence) in hypothesis testing and model specificity. It is hoped that the demonstrative analysis and discussion included in this study with the HVMs will also be useful for other ecological and evolutionary modeling of biodiversity and biogeography.

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

ZM designed and performed the study and wrote the manuscript.

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.