Introduction

An essential feature of microbial communities is their heterogeneous composition. A single environmental sample typically has a high richness, harboring hundreds to thousands of community members [1–3]. This high level of richness reaches an astronomical quantity at the global level, as scaling relationships and models of biodiversity predict upwards of one trillion (∼ 10¹²) species on Earth [4, 5]. Even experimental communities in laboratory settings with a single carbon source can harbor ≥ 40 community members, culminating in a total richness numbering in the hundreds among replicate communities (e.g., [6]). This richness contributes to the sheer diversity of microbial communities, a challenge for researchers attempting to identify the general principles that govern their dynamics and composition.

While richness estimates of microbial communities are undoubtedly high, the choice of assigning a community member to a given taxon remains intrinsically arbitrary. This arbitrariness remains regardless of whether the definition of a taxon is based on physiological attributes measured in the laboratory, entire genomes (i.e., metagenomics), or single-gene amplicon-based methods (i.e., 16S rRNA annotation). Despite their methodological differences, these approaches can all be viewed as different ways to cluster individuals within a community into groups. To contend with the sheer richness of microbial communities, researchers frequently rely on annotation-based approaches, i.e., by summing the abundances of species that belong to the same group at a given taxonomic scale (e.g., genus, family, etc.). This approach pares down communities to a size that is amenable for the visualization of individual groups and allows for questions of scale-dependent community reproducibility to be addressed [6–13].

The trend of performing analyses at a given taxonomic scale raises the question of how the composition of a community at one scale relates to that at another. To address these questions, researchers have examined the relationship between biodiversity measures at different scales in order to pare down the set of plausible ecological mechanisms that govern community composition. Specifically, recent efforts have found that microbial richness/diversity within a given taxonomic group (e.g., genus) is typically positively correlated with the richness/diversity among the remaining groups (e.g., family) [9, 14], an empirical pattern that aligns with the predictions of the Diversity Begets Diversity hypothesis (DBD) [15–17]. Evidence of the DBD hypothesis has historically been attributed to the construction of novel niches within a community through member interactions [15, 18], with similar mechanisms having been proposed to explain the existence of a positive relationship in microbial communities [14]. However, we still lack a quantitative understanding of how community composition at one scale should relate to that of another. Proceeding towards this goal requires two elements: 1) a systematic approach to grouping community members and 2) an appropriate null model for the composition of communities.

The operation of grouping the components of a system into a smaller number (e.g., merging read counts of OTUs to the family level in a community) is known in the physical sciences as coarse-graining. This formalism defines our systematic approach to grouping community members. While it is often not explicitly acknowledged as such, coarse-graining is a core concept in the microbial life sciences [19]. By smoothing over microscopic details at a lower level of biological organization in order to make progress at a higher level, the concept of coarse-graining has contributed towards the development of effective models of physiological growth [20, 21], evolutionary dynamics [22, 23], and the dependence of ecosystem properties on the diversity of underlying communities [24]. Coarse-graining has even been used to glean insight into the question of whether ”species” as a unit has meaning for microorganisms, as modeling efforts have found that the operation permits the delimitation of species when the resource preferences of community members are structured [25]. These theoretical and empirical efforts suggest that coarse-graining may provide an appropriate framework for investigating patterns of diversity and abundance within and between taxonomic scales of observation.

When evaluating the novelty of an empirical pattern it is useful to identify an appropriate null model for comparison [26–28]. Prior research efforts have demonstrated the novelty of the fine vs. coarse-grained relationship by contrasting inferences from empirical data with predictions obtained from the Unified Neutral Theory of Biodiversity (UNTB) [14, 29–33]. These predictions generally failed to reproduce slopes inferred from empirical data [14], implying that the fine vs. coarse-grained relationship represents a novel macroecological pattern that cannot be quantitatively explained by existing ecological models. However, the task of identifying an appropriate null model for comparison is not straightforward. Rather, the question of what constitutes an appropriate null model remains a persistent topic of discussion in community ecology [26, 34–37]. Here we take the view that a null model is appropriate for examining the relationship between two observables (e.g., community diversity at different scales) if it is capable of quantitatively predicting each observable (e.g., community diversity at one scale). By this standard, the UNTB is an unsuitable choice as a null as it generally fails to capture basic patterns of microbial diversity and abundance at any scale [38–40]. One relevant example is that the UNTB predicts that the distribution of mean abundances of species across sites is extremely narrow (i.e., converging to a delta distribution as the number of samples increases), whereas analyses suggest that empirical data tends to follow a broad lognormal distribution [40].

Contrastingly, recent efforts have determined that the predictions of a model of self-limiting growth with environmental noise, the Stochastic Logistic Model (SLM), is capable of quantitatively capturing multiple empirical macroecological patterns in microbial communities [40–44]. The stationary solution of this model predicts that the abundance of a given community member across sites follows a gamma distribution [40], a result that provides the foundation necessary to predict macroecological patterns among and between different taxonomic and phylogenetic scales.

In this study, we evaluated macroecological patterns of microbial communities across scales of evolutionary resolution. To limit potential biases that may result due to taxonomic annotation errors and to use all available data, we investigated the macroecological consequences of coarse-graining by developing a procedure that groups community members using the underlying phylogeny in addition to relying on taxonomic assignment. To ensure the generality of our findings and their commensurability with past research efforts, we used data from the Earth Microbiome Project (EMP), a public catalogue of microbial community barcode data. First, we assessed the extent that microbial diversity varies as the abundances of community members are coarse-grained by phylogenetic distance and taxonomic rank. The results of these analyses lead us to consider whether the predictive capacity of the gamma remains robust in the face of coarse-graining, a prediction that we quantitatively evaluated among community members and then extended to predict overall community richness and diversity. The accuracy of the gamma provided motivation to test whether the gamma was capable of predicting the relationship between fine and coarse-grained estimates, the empirical pattern that has been interpreted as evidence for the DBD hypothesis. Together, these analyses present evidence of the scale invariance of macroecological patterns in microbial communities as well as the applicability of the gamma distribution as a minimal model for evaluating macroecological patterns of microbial biodiversity.

Results

The macroecological consequences of phylogenetic and taxonomic coarse-graining

While microbial communities are often coarse-grained into higher taxonomic scales, their effect on measures of biodiversity and the underlying phylogeny are rarely examined. Before proceeding with the full analysis using public 16S rRNA amplicon data from the EMP, we elected to quantify the fraction of remaining community members across coarse-graining thresholds, a reflection of the extent that coarse-graining reduces global richness and the relation between taxonomic and phylogenetic coarse-graining. We first defined a coarse-grained group g as the set of OTUs that have the same assigned label in a given taxonomic rank out of G groups (e.g., Pseudomonas at the genus level) or are collapsed when the phylogeny is truncated by a given root-to-tip distance (Figs. 1, S1).The relative abundance of group g in site j is defined as .

Figure 1.

We found that coarse-graining had a drastic effect on the total number of community members within an environment, reducing global richness by ∼90% even at just the genus level (Fig. S2a). By coarse-graining over a range of phylogenetic distances, we saw that a fraction of coarse-grained community members comparable to that of genus-level coarse-graining occurs at a phylogenetic distance of ∼0.1 (Fig. S2b). This distance translates to only ∼3% of the total distance of the tree, meaning that the majority of OTUs are coarse-grained over a minority of the tree. This pattern is likely driven by the underlying structure of microbial phylogenetic trees, where most community members have short branch lengths [45]. This result suggests that while coarse-graining communities to the genus or family level substantially reduces global richness, it does so without coarse-graining the majority of the evolutionary history that is captured by the phylogeny. Assuming that phylogenies capture ecological changes that occur over evolutionary time, this detail implies that ecological divergence that is captured by the phylogeny should be retained even when communities are considerably coarse-grained.

With our coarse-graining procedures established, we proceeded with our macroecological investigation. Recent efforts have found that the distribution of abundances of a given ASV/OTU maintains a consistent statistically similar form across independent sites and time, a pattern known as the Abundance Fluctuation Distribution [40–42, 46, 47]. By coarse-graining empirical AFDs and rescaling them by their mean and variance (i.e., standard score), we saw that AFDs from the human gut microbiome retain their shape across phylogenetic scales (Fig. 2a). This pattern of invariance holds across environments for both phylogenetic and taxonomic coarse-graining (Figs. S3, S4), suggesting that AFDs can likely be described by a single probability distribution.

Figure 2.

It has been previously demonstrated that microbial AFDs are well-described by a gamma distribution that is parameterized by the mean relative abundance and the shape parameter (i.e., equal to the squared inverse of the coefficient of variation [40]). This distribution can be viewed as the stationary distribution of a Stochastic Logistic Model (SLM) of growth, a mathematical model that successfully captures macroecological patterns of microbial communities across both sites and time [40, 41, 44] (Eq. 5 in Materials and Methods).

Using this result, we determined whether the gamma sufficiently characterized coarse-grained AFDs. In order to accomplish this task, it is worth noting that we do not directly observe x_i. Rather, our ability to observe a community member is dependent on our sampling effort (i.e., total number of reads for a given site). To account for sampling, one can derive a form of the gamma distribution that explicitly accounts for the sampling process, obtaining the probability of obtaining n reads out of N total reads belonging to a community member Pr (Materials and Methods, [40]). Given that n = 0 for a species we do not observe, we can define the fraction of M sites where a species was observed (i.e., occupancy, o_i) as

We can then compare this prediction to observed estimates of occupancy to assess the accuracy of the gamma distribution across coarse-grained thresholds. We found that Eq. 1 generally succeeded in predicting observed occupancy across phylogenetic and taxonomic scales for all environments (Figs. S5, S6). We then determined whether the gamma is capable of predicting the relationship between macroecological quantities. One such relationship is that the occupancy of a species should increase with its mean abundance, known as the abundance-occupancy relationship [48]. This pattern has been found across microbial systems [49–51] and can be quantitatively predicted using the gamma [40]. We see that this relationship is broadly captured across taxonomic and phylogenetic scales for all environments (Figs. 2b, S7 S8). This result implies that the ability to observe a given taxonomic group is primarily determined by its mean abundance across sites and the sampling effort within a site [40].

To quantitatively assess the gamma we calculated the relative error of our occupancy predictions (Eq. 10) for all coarse-graining thresholds. We found that the mean logarithm of the error only slightly increased for the initial taxonomic and phylogenetic scales, where it then exhibited a sharp decrease across environments (Fig. S9a,b). The error then only began to decrease once the community became highly coarse-grained, harboring a global richness (union of all community members in all sites for a given environment) < 20. This result means that, if anything, the accuracy of the gamma distribution only improves with coarse-graining.

Reconciling coarse-graining and the predictions of the gamma distribution

The consistent predictive success of the gamma distribution in the face of coarse-graining raises the question of why it remains a sufficient null model. The sum of independent gamma distributed random variables only returns a gamma if all variables have identical rate parameters , a requirement that microbial communities clearly do not meet since they typically harbor broad mean abundance distributions. Given that a gamma AFD cannot predict the distribution of correlations between AFDs [40], it is first worth examining whether the degree of dependence between AFDs shapes coarse-grained variables. We first consider the relation between the variance of the sum and the sum of variances.

By plotting against across coarse-grained thresholds, we found that the contribution of covariance was weak, meaning that the statistical moments at higher scales can be approximated by those at lower scales (Figs. S10, S11). This result is consistent with previous efforts demonstrating that the strongest correlations between AFDs are typically concentrated among pairs of closely related community members (i.e., low phylogenetic distance) [52], implying that the effect of correlation should dissipate when communities are coarse-grained by their taxonomy or phylogeny. Given that the variance of the sum can be approximated by the sum of the variances and that, by definition, the mean of a sum is the sum of the means, it is reasonable to propose that the statistical moments of coarse-grained AFDs are sufficient to characterize the distribution.

Finally, while there is no general closed-form solution for the sum of independent gamma distributed random variables with different rate parameters (equivalent to considering the convolution of many AFDs with different carrying capacities), continued progress has been made towards obtaining suitable approximations [53–57]. This body of work includes an analysis demonstrating that a single gamma distribution can provide a suitable approximation to the distribution of the sum of many gamma random variables with different rate parameters [58]. In summary, the gamma distribution successfully captures patterns of biodiversity under taxonomic and phylogenetic coarse-graining because the correlation between community members is typically low and the sum of multiple gamma distributions can be approximated by a single gamma distribution.

Predicting measures of richness and diversity within a coarse-grained scale

Given that the presence or absence of a community member is used to estimate community richness, a measure previously used to make claims about patterns of microbial diversity across taxonomic scales [14], we can visualize the sufficiency of the gamma by predicting the mean richness within an environment at a given coarse-grained threshold (Eq. 11). Likewise, we can use the entirety of the distribution of read counts to predict the diversity within a site, a measure that reflects richness as well as the distribution of abundances within a community (Eq. 12), analytic predictions that we validated through simulations (Fig. S12).

Focusing on the human gut microbiome as an example, we found that we can predict the typical richness of a community across phylogenetic scales using the gamma distribution (Fig. 3a). This success is repeated when we examined the predicted diversity (Fig. 3b). By examining all nine environments, we see that despite the dissimilarity in environments we were able to predict mean richness and diversity in the face of coarse-graining (Fig. 3c,d). The results of this analysis suggest that the composition of microbial communities remains largely invariant in the face of coarse-graining and that the gamma distribution remains a suitable null model for predicting average community measures across coarse-grained scales. Identical results were obtained for taxonomic coarse-graining (Fig. S13).

Figure 3.

Turning to higher-order moments, we examined the variance of richness and diversity across sites. Using a similar approach that was applied to the mean, we derived analytic predictions for the variance (Eq. 15). With the human gut as an example, we see that analytic predictions typically fail to capture estimates of variance obtained from empirical data for phylogenetic coarse-graining (Fig. 4a,b). This lack of predictive success was consistent across environments (Fig. 4c,d), implying that a model of independent species with gamma distributed abundances is insufficient to capture the variance of measures of biodiversity. A major assumption made in our derivation was that species are independent, an assumption that is almost certainly wrong given that the gamma distribution has been previously demonstrated to be unable to capture the empirical distribution of correlations in the AFDs of community members [40]. To attempt to remedy this failed prediction, we again turned to the law of total variance by estimating the covariance of richness and diversity from empirical data and adding the covariance to the predicted variance for each measure. We found that the addition of this empirical estimate was sufficient to predict the observed variance in the human gut (Fig. 4a,b) as well as across environments (Fig. 4e,f), implying that the underlying model is fundamentally correct for predicting the first moment of measures of biodiversity but cannot capture the correlations necessary to explain higher statistical moments such as the variance. Identical results were again obtained with taxonomic coarse-graining (Fig. S14).

Figure 4.

Predicting patterns of richness and diversity between coarse-grained scales

Our predictions of the statistical moments of richness and diversity using the gamma distribution provided the foundation necessary to investigate macroecological patterns between different coarse-grained scales. One such prominent pattern is the relationship between the fine-grain richness/diversity within a given coarse-grained group vs. the coarse-grained richness/diversity among all remaining groups (e.g., the number of classes within Firmicutes vs. the number of phyla excluding the phylum Firmicutes), a pattern that has been purported to reflect the existence of DBD in microbial systems. Before continuing, we note that the acronym DBD technically refers to the hypothesis that such positive relationships reflects the existence of ecological interactions through which coarse-grained diversity bolsters the accumulation of fine-grained diversity (e.g., niche construction [16, 59]). Since we are primarily interested in the predictive power of an empirically-validated null model of biodiversity, we distinguish between DBD as a hypothesis and DBD as an empirical pattern by referring to the slope as the fine vs. coarse-grained relationship throughout the remainder of this manuscript.

The fine vs. coarse-grained relationship can be quantified as the slope of the relationship between the fine-grained richness within a given coarse-grained group g (S_g,m) and the richness in the remaining G − 1 coarse-grained groups: S_g,m ∝ αS_G\g,m, where G \ g denotes the exclusion of group g and α is the slope of the relationship. This formulation was proposed in Madi et al. and to ensure commensurability we adopted it here [14]. Using Eq. 1, we can define each of these estimators in terms of the sampling form of the gamma while accounting for sampling

Similarly, we can use Eq. 12 to derive predictions for fine and coarse-grained diversity.

By repeating this calculation for all M sites, we obtained vectors of coarse and fine-grained richness estimates for group g from which we can infer the fine vs. coarse-grained slope through ordinary least squares regression. By repeating this process for all G groups we obtained a distribution of slopes that can be directly compared to those obtained from empirical data.

Before performing a direct comparison, we first note the features of the empirical slopes and how they pertain to the predictions we obtained. By examining the distribution of empirical slopes pooled over all coarse-graining thresholds for each environment, we see that they are rarely less than zero (Figs. 5a, S19a). The few negative slopes inferred from empirical data were extremely small, having absolute values < 10⁻⁴ and can be treated as zeros. Furthermore, the distribution of slopes follows the same form across environments, suggesting that the fine vs. coarse-grained slope reflects a general feature of community sequence data rather than the ecology of specific environments. Like the empirical slopes, the gamma distribution virtually always predicted a positive slope for all environments for both taxonomic and phylogenetic coarse-graining. This paucity of negative slopes suggests that the prediction of the alternative to the DBD hypothesis, the Ecological Controls hypothesis [60], is virtually absent in empirical data and cannot be generated from an empirically validated null model of microbial biodiversity.

Figure 5.

However, only observing positive slopes does not necessarily provide support to the DBD hypothesis. A direct comparison of slopes predicted from the gamma distribution to those inferred from empirical data is necessary to determine whether the predictions of DBD lie outside what can be reasonably captured by an interaction-free model such as the SLM. To evaluate the novelty of the fine vs. coarse-grained slopes we compared the values of observed slopes to those obtained from the interaction-free SLM. We found that the predictions of the gamma distribution closely matched the observed slopes across environments for both taxonomic and phylogenetic coarse-graining (Figs. S15, S16). We can consolidate these results by taking the average slope for a given coarse-grained level, from which we see that the mean slope predicted by the gamma distribution does a reasonable job capturing empirical slopes across environments (Figs. 5b, S19b). These results indicate that we should expect to see a positive relationship between richness estimates at different scales and that the relationships we observe can be quantitatively captured by a gamma distributed AFD. It is worth noting that the fine vs. coarse-grained slope could be sufficiently predicted even though the gamma distribution only succeeded at predicting average richness, suggesting that higher order statistical moments, and by extension interactions between community members, are unnecessary to quantitatively capture the positive relationship observed between fine and coarse-grained estimates of richness.

While richness is a widespread and versatile estimator that is commonly used in community ecology, it can neglect a considerable amount of information by focusing on presences and absences instead of the entirety of the distribution of abundances. To rigorously test the predictive power of the gamma distribution it is necessary to evaluate the fine vs. coarse-grained relationship for diversity. We again found that disparate environments had similar distributions of slopes from empirical data (Figs. 5c, S19c), suggesting that the slope of the relationship is likely a general property of microbial communities rather than an environment-specific pattern. However, unlike richness, diversity predictions obtained from the gamma distribution generally failed to capture observed slopes, as the squared correlation between observed and predicted slopes can be less than that of richness by over an order of magnitude (Figs. 5d, S17, S18, S19d). Here we see where the predictions of an interaction-free SLM succeeded and failed to predict observed macroecological patterns.

Given that the gamma distribution failed to predict the observed diversity slope, it is worth evaluating whether additional features could be incorporated to generate successful predictions. A notable omission is that there is an absence of interactions between community members in the SLM, meaning that we were unable to predict correlations between community member abundances. However, while considerable progress has been made (e.g., [11]), predicting the observed distribution of correlation coefficients between community members while accounting for sampling remains a non-trivial task. Given that the gamma distribution succeeded at predicting other macroecological patterns, we elected to perform a simulation where a collection of sites was modeled as an ensemble of communities with correlated gamma-distributed AFDs with the means, variances, correlations, and total depth of sampling set by estimates from empirical data (Materials and methods). By including correlations between AFDs into the simulations, the statistical outcome of ecological interactions between community members, we were able to largely capture the observed fine vs. coarse-grained diversity slopes (Figs. 6, S20, S21, S22). These results suggest that rather than diversity at a fine-scale begetting diversity at a coarse-scale, the correlations that exist at a fine-scale do not dissipate under coarse-graining, resulting in a positive relationship between measures of diversity at different scales.

Figure 6.

Discussion

The results of this study demonstrate that macroecological patterns in microbial communities remain largely invariant across taxonomic and phylogenetic scales. By focusing on the predictions of the SLM, an interaction-free model of microbial growth under environmental fluctuations, we were able to evaluate the extent that measures of biodiversity can be predicted under coarse-graining. We were largely able to predict said measures using the same model with zero free parameters (i.e., no statistical fitting) across scales, implying that certain macroecological patterns of microbial communities remain self-similar across taxonomic and phylogenetic scales. Building off of this result, we investigated the dependence of community measures between different degrees of coarse-graining, a pattern that has been formalized as the Diversity Begets Diversity hypothesis [14, 15]. The prediction derived from the sampling form of the gamma distribution quantitatively captured the observed slopes of the fine vs. coarse-grained relationship for richness, while it failed to capture the slope of diversity. We showed that introducing correlations between abundance fluctuations in the SLM allow to recover the prediction for the slopes of diversity.

Our richness results complement past work which demonstrated that occupancy, the constituent of richness, is highly dependent on sampling depth (i.e., total read count) and the mean abundance of a community member [40]. This past work and the relationships between the mean abundance and occupancy evaluated in this manuscript demonstrate that occupancy alone is unlikely to contain ecological information that is not already captured by the distribution of abundances across hosts (i.e., the AFD). Our analyses of the relationship between fine and coarse-grained richness support this conclusion, as predictions derived from a gamma distribution quantitatively captured the observed slope. The success of an interaction-free model in predicting fine vs. coarse-grained slopes is an indictment of the appropriateness of estimators that rely solely on the presence of a community member for identifying novel macroecological patterns, a measure that has been used to bolster support for the DBD hypothesis at the level of 16S rRNA amplicons as well as strains [14, 61]. Rather, estimates of richness harbor little information about the dynamics of a community across taxonomic and phylogenetic scales that is not already captured by the sampling form of the gamma. Contrasting with richness, the predictions of diversity from the gamma sampling distribution were unable to capture fine vs. coarse-grained relationships from empirical data. Given that measures of diversity incorporate information about the richness and evenness of a community [62], the comparative deficiency of our predictions for fine vs. coarse-grained diversity suggests that forms of the SLM that neglect interactions between community members cannot capture coarse-graining relationships that depend on the evenness of the distribution of abundances.

Macroecological patterns are not imbued with mechanistic explanation [63]. Rather, the onus is on the investigator to identify plausible mechanisms. Often in ecology this task is made easier by eliminating mechanisms that are incapable of producing the observed pattern, that is, identifying an appropriate null. The novelty of the fine vs. coarse-grained relationship was previously assessed using a null model which assumed that community dynamics were primarily driven by demographic noise (i.e., the UNTB) [14, 32]. Empirical patterns of microbial abundance cannot be reasonably captured by such models, making predictions obtained from the UNTB invalid for evaluating the novelty of microbial macroecological patterns. In contrast, models that combine self-limiting growth with environmental noise reproduce several empirical patterns, making the SLM an appropriate choice for evaluating the novelty of fine vs. coarse-grained relationships [40, 44]. This is not a trivial detail, as there is historical precedence on the need to identify an appropriate null in order to investigate how fine and coarse-grained measures of biodiversity relate to one another, as one of the earliest adoptions of null model analysis in ecology was done to investigate the ratio of species to genera in a community [64, 65].

In this study, the predictions of the sampling form of the gamma considerably improved when correlations between community members were included. This result suggests that rather than exclusively pointing to niche construction as previously suggested [14], any ecological mechanism that can capture the observed distribution of correlation coefficients is a plausible candidate. Given that models of consumer-resource dynamics have succeeded in capturing macroecological patterns [66, 67], including quantitatively predicting the distribution of correlation coefficients [11], it is reasonable to suggest that such mechanisms are ultimately responsible for the relationship between fine and coarse-grained measures of diversity and can be reduced to phenomenological models such as the SLM. Indeed, experimental investigations of the slope evaluated here have found that positive slopes exist in artificial communities maintained in a laboratory setting, where the strength of the correlation between fine and coarse-grained scales is driven by the secretion of secondary metabolites [9]. This phenomenon, known as cross-feeding, can be viewed as a mechanism that is compatible with the concept of niche construction [16] and the original interpretation of Madi et al. [14].

In the interest of providing macroecological insight into the DBD hypothesis, we solely focused on coarse-graining procedures that relied on phylogenetic reconstruction and taxonomic assignment. However, it is worth noting that it is also possible to coarse-grain community members by the strength of their correlations (i.e., sum the abundances of each pair of community members with the strongest correlation in AFDs). This procedure has been named the phenomenological renormalization group method due to its ability to identify if and where a system is stable despite knowing little about the system’s dynamics (i.e., fixed points in nonlinear systems) [68, 69].

However, given that the AFD correlation between two community members is often inversely related to their phylogenetic distance, such an analysis would likely be redundant, as coarse-graining based on the strength of correlation would effectively coarse-grain the most closely related community members [52].

A major goal of this study was to evaluate the novelty of macroecological patterns that were used to bolster support for the DBD hypothesis. We used the same dataset in order to ensure generality and commensurability with past research efforts. However, it is worth inspecting how the use of a global survey dataset constrains the inferences one can make. Throughout this study we implicitly assumed that an ensemble approach is valid, meaning that we viewed different sites/hosts as virtual copies of a given environment. This assumption is valid for time-series studies where the distribution of microbial abundances remains stationary with respect to time [70], as the stationary solution of the SLM has successfully characterized microbial community time-series at both the level of OTUs [40] and strains [47]. So, in communities that are stationary with respect to time, our fine vs. coarse-grained relationship results will remain valid.

Materials and methods

Data acquisition and processing

To ensure that our analyses were generalizable across ecosystems and commensurate with prior DBD investigations, we used amplicon sequence data from the V4 region of the 16S rRNA gene generated and curated by the Earth Microbiome Project [1, 14]. We restricted our analysis to the quality control (QC)-filtered subset of the EMP, which was annotated using the closed-reference database SILVA [71] and consists of 96 studies culminating in 23,828 total samples with each processed sample having ≥10,000 reads. We downloaded the public Silva reference tree for OTUs with 97% similarity 97 otus.tre from the EMP database. We identified nine heavily sampled environments in the metadata file emp qiime mapping qc filtered.tsv and selected 100 random samples from each environment.

We briefly note that our occupancy and richness predictions depend on the form of the gamma distribution that explicitly accounts for sampling as a multinomial process. The multinomial distribution describes the probability of sampling n reads given a relative abundance of x and total read count N with replacement, a process we can model as the Poisson limit of a binomial sampling process for individual community members. Given this choice and the past success of the gamma distribution, we deviated from past analyses by electing to not sub-sample read counts to the same depth across samples, as the process of sampling without replacement would bias the sampling distribution for rare community members [14].

Coarse-graining protocol

Taxonomic coarse-graining was performed as the summation of the abundances of all OTUs within a given taxonomic group. We removed taxa with indeterminate labels to prevent potential biases due to taxonomic misassignment, (e.g., ”uncultured”, ”ambiguous taxa”, ”candidatus”, ”unclassified”, etc.). Manual inspection of EMP taxonomic annotations revealed a low number of OTUs that had been assigned the taxonomic label of their host (e.g., Arachis hypogaea (peanut)). These marked OTUs were removed from all downstream analyses.

Phylogenetic coarse-graining was performed using the phylogenetic tree provided by SILVA 123 97 otus.tre in the EMP release. Each internal node of a phylogenetic tree was collapsed if the mean branch lengths of its descendants was less than a given distance. All phylogenetic operations were performed using the Python package ETE3 [72].

Deriving biodiversity measure predictions

While the gamma distribution as the stationary solution of the SLM and the sampling form of the gamma have been previously derived [40], we briefly outline relevant derivations here for the convenience of the reader before deriving the predicted richness and diversity of a community. We define the SLM as the following Langevin equation

Here τ_i, K_i, and σ_τi represent the timescale of growth, the carrying capacity, and the coefficient of variation of growth rate fluctuations, respectively. Multiplicative environmental noise is captured by the product of a linear frequency term, the coefficient of variation of growth rate fluctuations, and a Brownian noise term η(t) that introduces stochasticity into the equation. The expected value of η(t) is ⟨η(t)⟩= 0 [73].

The dependence of η(t^′) at time t^′ on an earlier time η(t) is defined as ⟨η(t)η(t^′)⟩ = δ(t −t^′) [73]. This standard definition means that if the noise term is shifted in time, it has zero correlation with itself. We briefly note that because DBD patterns were originally investigated in Madi et al. using an ensemble of sites that belong to the same type of ecosystem rather than the timeseries of a single site [14], the gamma distribution alone does not prove the validity of the SLM nor does it prove alternatively formulated stochastic differential equations of ecology that also predict a gamma distribution (e.g., [74]). However, given that the SLM has successfully characterized the temporal dynamics of microbial communities, we believe that this model is an appropriate formulation for investigating DBD patterns [40, 44, 47].

In contrast to the SLM, macroecological predictions can be derived from the UNTB. There are many forms of the UNTB, but the novelty of observed fine vs. coarse-grained relationships was assessed using a form of the UNTB that predicts that the distribution of community member abundances within a given site follows a zero-sum multinomial distribution [14, 32]. Formulations of the UNTB generally seek to capture ecological dynamics in a community without fixed carrying capacities that are driven by demographic noise, so for the purpose of this study the SLM can be viewed in contrast to the following Langevin model of neutral ecology: .

The stationary distribution of the SLM can be derived using the Itô ↔ Fokker-Planck equivalence and taking the stationary solution [40, 75], resulting in the gamma distributed AFD. Through the SLM, we can define the mean relative abundance and its squared inverse coefficient of variation as and , respectively. These are parameters that are estimated from the empirical data and will be used below to obtain predictions. Using these definitions and the stationary distribution of Eq. 5, we obtain the gamma distribution

When we sequence microbial communities, we do not directly observe abundances. Rather, we obtain read counts from sequencing. Therefore, it is necessary to account for the reality of sampling when we apply Eq. 6 to empirical data. We can account for sampling by first assuming that the probability of observing a single community member can be modeled as a binomial sampling process. Given that the total number of reads is typically large (N ≫ 1) and the typical relative abundance of a community member is much smaller than one (x_i ≪1), the binomial can be approximated as a Poisson sampling process with the following probability of sampling n reads

This formulation of the sampling process is convenient, as it can be used to obtain an analytic solution for the probability of observing n reads given and β_i, the parameters we estimate from the data. This distribution can be obtained solving the convolution of the Poisson and the gamma [40]. The result allows us to calculate the probability of obtaining n reads out of a total sampling depth of N for the ith OTU as

We will be using this equation to obtain predictions of measures of biodiversity. First, noticing that the probability of a community member’s absence is the complement of its presence, we can define the expected occupancy of a community member across M sites as

The success of our predictions was assessed using the relative error.

Using the definition of occupancy from the sampling form of the gamma, we derive the expected richness of a community as

Similarly, we can derive the expected diversity as

In physics parlance, these predictions neglect interactions between community members, also known as mean-field predictions. We can calculate the mean-field prediction of Eq. 11 from empirical data. However, there is no analytic solution for the integral inside the sum of Eq. 12. To calculate ⟨H⟩, we performed numerical integration on each integral for each taxon in each sample at a given coarse-grained resolution using the quad() function from SciPy.

To predict the variance of each measure we derived the expected value of the second moment, assuming independence among community members. We derive the second moments of richness and diversity

where δ_i,j is the Kronecker delta.

By performing an analogous series of operations, we obtain the expected value of the second moment for diversity.

Where the expected value of the second moment of the diversity term is defined as . From which we obtain the expected value of the variance

Fine vs. coarse-grained slope inference

In order to predict the relationship between the measures within a coarse-grained group and that among all remaining groups, we calculated a vector of predicted richness or diversity estimates for all sites using Eq. 3a or Eq. 4a within a given coarse-grained group and 3b or Eq. 4b among the remaining groups. This “leave-one-out” procedure was originally implemented in Madi et al., where the authors examined the slope of fine vs. coarse-grained measures of diversity as a sliding window across taxonomic ranks with both the fine and coarse scales increasing with each rank (e.g., genus:family, family: order, etc.) [14]. To maintain consistency, we used the same definition for our predictions. We also extended the definition to the case of phylogenetic coarse-graining, where we compared fine and coarse scales using different phylogenetic distances while retaining the same ratio (e.g., 0.1:0.3, 0.3:0.5, etc.).

Slopes were estimated using ordinary least squares regression with SciPy. Throughout the manuscript the success of a prediction was evaluated by calculating its relative error as follows. We only inferred the slope if a fine-grained group had at least five members. We only examined the slopes of a given coarse-grained threshold if at least three slopes could be inferred.

Simulating communities of correlated gamma-distributed AFDs

Correlated gamma-distributed AFDs were simulated by performing inverse transform sampling. For each environment with M sites, an M × S_obs matrix Z was generated from the standard Gaussian distribution using the empirical S_obs × S_obs correlation matrix of empirical relative abundances. The cumulative distribution U = Φ (Z)_Gaus. was obtained and a matrix of SADs was obtained using the point percentile function of the gamma and empirical distribution of mean relative abundances and the squared inverse coefficient of variation of abundances: , β = β₁, β₂, ⋯, . To simulate the process of sampling, each SAD of the resulting M × S_obs matrix of true relative abundances was sampled using a multinomial distribution with the empirical distribution of total read counts.

Data and code availability

All sequencing data used in this study was obtained from the EMP (URL: ftp://ftp.microbio.me/emp/release1/). Processed data are available on Zenodo, DOI: 10.5281/zenodo.7692046. All code written for this study is available on GitHub under a GNU General Public License: macroeco phylo

Acknowledgements

This work was supported by the NSF Postdoctoral Research Fellowships in Biology Program under Grant No. 2010885 (W.R.S.).

Author contributions

W.R.S. and J.G. conceptualized the project, completed the derivations, and wrote the manuscript. W.R.S. performed all analyses.