Investigating macroecological patterns in coarse-grained microbial communities using the stochastic logistic model of growth

  1. William R Shoemaker  Is a corresponding author
  2. Jacopo Grilli
  1. Quantitative Life Sciences, The Abdus Salam International Centre for Theoretical Physics (ICTP), Italy


The structure and diversity of microbial communities are intrinsically hierarchical due to the shared evolutionary history of their constituents. This history is typically captured through taxonomic assignment and phylogenetic reconstruction, sources of information that are frequently used to group microbes into higher levels of organization in experimental and natural communities. Connecting community diversity to the joint ecological dynamics of the abundances of these groups is a central problem of community ecology. However, how microbial diversity depends on the scale of observation at which groups are defined has never been systematically examined. Here, we used a macroecological approach to quantitatively characterize the structure and diversity of microbial communities among disparate environments across taxonomic and phylogenetic scales. We found that measures of biodiversity at a given scale can be consistently predicted using a minimal model of ecology, the Stochastic Logistic Model of growth (SLM). This result suggests that the SLM is a more appropriate null-model for microbial biodiversity than alternatives such as the Unified Neutral Theory of Biodiversity. Extending these within-scale results, we examined the relationship between measures of biodiversity calculated at different scales (e.g. genus vs. family), an empirical pattern previously evaluated in the context of the Diversity Begets Diversity (DBD) hypothesis (Madi et al., 2020). We found that the relationship between richness estimates at different scales can be quantitatively predicted assuming independence among community members, demonstrating that the DBD can be sufficiently explained using the SLM as a null model of ecology. Contrastingly, only by including correlations between the abundances of community members (e.g. as the consequence of interactions) can we predict the relationship between estimates of diversity at different scales. The results of this study characterize novel microbial patterns across scales of organization and establish a sharp demarcation between recently proposed macroecological patterns that are not and are affected by ecological interactions.

eLife assessment

This valuable study considers empirical macroecological patterns in microbiome data across multiple taxonomic scales. The work convincingly shows that the Stochastic Logistic Growth model is a more appropriate choice of null model than the neutral theory of biodiversity. The work will be of particular interest to microbial ecologists.


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 (Thompson et al., 2017; Shoemaker et al., 2017; Barberán et al., 2014). 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 (∼1012) species on Earth (Locey and Lennon, 2016; Lennon and Locey, 2020). 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. Dal Bello et al., 2021). 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 exists 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 community members 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 (Louca et al., 2016; Goldford et al., 2018; Estrela et al., 2022; Estrela et al., 2021; Dal Bello et al., 2021; Ho et al., 2022; Good and Rosenfeld, 2022; Tian et al., 2020).

This movement towards performing analyses of diversity at various taxonomic scales 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 ecological mechanisms that plausibly 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) (Madi et al., 2020; Estrela et al., 2022), an empirical pattern that aligns with the predictions of the Diversity Begets Diversity (DBD) hypothesis (Whittaker, 1972; San Roman et al., 2018; Maynard et al., 2017). Evidence of the DBD hypothesis has historically been attributed to the construction of novel niches within a community through member interactions (Calcagno et al., 2017; Whittaker, 1972), with similar mechanisms having been proposed to explain the existence of a positive relationship in microbial communities (Madi et al., 2020). 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 (Good and Hallatschek, 2018). 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 (Scott et al., 2014; Jun et al., 2018), evolutionary dynamics (Schweinsberg, 2003; Desai et al., 2013), and the dependence of ecosystem properties on the diversity of underlying communities (Moran and Tikhonov, 2022). 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 distinct taxonomic groups when the resource preferences of community members are structured (Tikhonov, 2017). 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 (O’Dwyer et al., 2017; McGill, 2010; Harte, 2011). 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) (Hubbell, 2011; Volkov et al., 2003; Azaele et al., 2016; Madi et al., 2020; Alonso and McKane, 2004; Azaele et al., 2006). These predictions generally failed to reproduce slopes inferred from empirical data (Madi et al., 2020), implying that the fine vs. coarse-grained relationship represents a novel macroecological pattern that cannot be quantitatively explained by existing null models of ecology. 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 (Simberloff, 1983; Harvey et al., 1983; Gotelli and Graves, 1996; Gotelli and Ulrich, 2012; O’Dwyer et al., 2017). 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 was 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 (Li and Ma, 2016; Harris et al., 2017; Grilli, 2020). One relevant example is that the UNTB predicts that the distribution of mean abundances of community members across sites is extremely narrow (i.e. converging to a delta distribution as the number of sites increases), whereas empirical data tends to follow a broad lognormal distribution (Grilli, 2020). Contrastingly, recent efforts have determined that the predictions of a model of self-limiting growth with environmental noise, the SLM, is capable of quantitatively capturing multiple empirical macroecological patterns in observational and experimental microbial communities (Grilli, 2020; Zaoli and Grilli, 2021; Zaoli et al., 2022; Descheemaeker et al., 2021; Descheemaeker and de Buyl, 2020; Shoemaker et al., 2023c; Lim et al., 2023). The stationary solution of this model predicts that the abundance of a given community member across sites follows a gamma distribution (Grilli, 2020), 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. We used data from the Earth Microbiome Project (EMP), a public catalog of microbial community barcode data, to ensure the generality of our findings and their commensurability with past research efforts. 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 distribution remained robust under 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 distribution provided the necessary motivation to test whether the gamma distribution 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, the stationary distribution of the SLM, as a null model for evaluating the novelty of macroecological patterns of microbial biodiversity.


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 the 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 (Figure 1, Figure 1—figure supplement 1). The relative abundance of group g in site j is defined as xg,j=igxi,j.

Figure 1 with 2 supplements see all
The process of coarse-graining abundances using the phylogeny.

Taxonomic assignment in 16S rRNA amplicon sequence data provided the opportunity to investigate how properties of communities vary at different taxonomic scales. The most straightforward means of coarse-graining here is to sum the abundances of OTUs/ASVs that belong to the same taxonomic group. Amplicon data-based studies provide information about the shared evolutionary history of community constituents, information that can be leveraged by the construction of phylogenetic trees. A coarse-graining procedure can be defined that is analogous to one based on taxonomy, where a phylogenetic root-to-tip distance is chosen and terminal nodes are collapsed if their distance to a common ancestor is less than the prescribed distance.

We found that even minor degrees of 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 (Figure 1—figure supplement 2a). By coarse-graining over a range of phylogenetic distances, we found that the fraction of coarse-grained community members comparable to that of genus-level coarse-graining occurred at a root-to-tip phylogenetic distance of ∼0.1 (Figure 1—figure supplement 2b). This distance translated to only ∼3% of the total distance of the tree, meaning that the majority of OTUs were coarse-grained over a minority of the tree. This pattern was likely driven by the underlying structure of microbial phylogenetic trees, where most community members have short branch lengths (O’Dwyer et al., 2015). 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 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 maintained a consistent statistically similar form across independent sites and time, a pattern known as the Abundance Fluctuation Distribution (Grilli, 2020; Zaoli and Grilli, 2021; Zaoli et al., 2022; Shoemaker, 2023a; Wolff et al., 2023). By coarse-graining empirical AFDs and rescaling them by their mean and variance across sites (i.e. standard score), we found that AFDs from the human gut microbiome retained their shape across phylogenetic scales (Figure 2a). This pattern of invariance held across environments for both phylogenetic and taxonomic coarse-graining (Figure 2—figure supplement 1, Figure 2—figure supplement 2), suggesting that empirical AFDs can likely be described by a single probability distribution.

Figure 2 with 13 supplements see all
The shape of the AFD remained qualitatively invariant under coarse-graining.

(a) Under phylogenetic coarse-graining the general shape of the AFD for OTUs that were present in all sites (i.e. an occupancy of one) remained qualitatively invariant. (b) Similarly, the shape of the relationship between the mean coarse-grained abundance across hosts and occupancy across sites did not tend to vary. Predictions obtained from the gamma distribution are capable of capturing the relationship between the mean abundance and occupancy, suggesting that the gamma distribution remains a useful quantitative null model under coarse-graining. All data in this plot is from the human gut microbiome.

It has been previously demonstrated that empirical microbial AFDs are well-described by a gamma distribution that is parameterized by the mean relative abundance x¯i and the shape parameter βi=x¯i2σi2 (equal to the squared inverse of the coefficient of variation Grilli, 2020). This distribution can be viewed as the stationary distribution of a SLM of growth, a mathematical model that successfully captures macroecological patterns of microbial communities across both sites and time (Grilli, 2020; Descheemaeker and de Buyl, 2020; Zaoli and Grilli, 2021; Equation 5 in Materials and methods).

Using this result, we determined whether the gamma distribution sufficiently characterized coarse-grained AFDs. In order to accomplish this task, it is worth noting that we do not directly observe xi. Rather, our ability to observe a community member is dependent on 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 (Materials and methods, Grilli, 2020). Given that n=0 for a community member, we do not observe, we defined the fraction of M sites where a community member was observed (i.e. occupancy, oi) as

(1) oi=11Mm=1MP(0|Nm,x¯i,βi)

We then compared this prediction to observed estimates of occupancy to assess the accuracy of the gamma distribution across coarse-grained thresholds. We found that Equation 1 generally succeeded in predicting observed occupancy across phylogenetic and taxonomic scales for all environments (Figure 2—figure supplement 3, Figure 2—figure supplement 4). We then determined whether the gamma distribution was capable of predicting the relationship between macroecological quantities. One such relationship is that the occupancy of a community member should increase with its mean abundance, known as the abundance-occupancy relationship (Gaston et al., 2000). This pattern has been found across microbial systems (Shade et al., 2018; Sloan et al., 2007; Burns et al., 2016) and can be quantitatively predicted using the gamma distribution (Grilli, 2020). We see that this relationship is broadly captured across taxonomic and phylogenetic scales for all environments (Figure 2b, Figure 2—figure supplement 5, Figure 2—figure supplement 6). This result implies that the ability to observe a given taxonomic group was primarily determined by its mean abundance across sites and the sampling effort within a site, regardless of one’s scale of observation. In contrast, under the assumption of demographic indistinguishability under the UNTB we would expect the mean abundance distribution to be extremely narrow, following a delta distribution. Under the SLM, the variation in mean relative abundances we observed implies that the carrying capacities of community members vary over multiple orders of magnitude. We also note that at high mean abundances our predictions show slight variation, which is likely driven by variation in the shape parameter β (Figure 2b). In contrast with these results, the gamma distribution was unable to predict the variance of occupancy under both taxonomic and phylogenetic coarse-graining (Equation 11; Figure 2—figure supplement 7, Figure 2—figure supplement 8), the implications of which we will address in a later section.

To quantitatively assess the accuracy of the gamma distribution we calculated the relative error of our mean occupancy predictions (Equation 12) 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 (Figure 2—figure supplement 9). 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 improved with coarse-graining.

Reconciling coarse-graining and the predictions of the gamma distribution

The consistent predictive success of the gamma distribution under 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 through analytic calculation if all random variables have identical rate parameters (βi/x¯i=β/x¯), 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 (Grilli, 2020), 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.

(2) Var(iStotxi)=iStotVar(xi)+2i<jCov(xi,xj)

By plotting Var(iStotxi) against iStotVar(xi) across coarse-grained thresholds, we found that the contribution of covariance to individual coarse-grained taxa was weak, suggesting that the statistical moments at higher scales can be approximated by those at lower scales (Figure 2—figure supplement 10, Figure 2—figure supplement 11). Similar conclusions can be drawn by plotting the variances as a ratio, with slight deviations above a ratio of one, suggesting that coarse-grained variance was slightly higher (Figure 2—figure supplement 12, Figure 2—figure supplement 13). These results are 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) (Sireci et al., 2023), implying that the effects of correlation should dissipate when communities are coarse-grained. 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 we know of 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), progress has been made towards obtaining suitable approximations (Stewart et al., 2007; Murakami, 2015; Hu et al., 2020; Behme and Bondesson, 2017; Barnabani, 2017). 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 (Covo and Elalouf, 2014). In summary, the gamma distribution appears to successfully captures patterns of biodiversity under taxonomic and phylogenetic coarse-graining because 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 (Madi et al., 2020), we can visualize the sufficiency of the gamma distribution by predicting the mean richness within an environment at a given coarse-grained scale (Equation 13a). 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 (Equation 14a), analytic predictions that we validated through simulations (Figure 3—figure supplement 1). We note that we observe consistent deviations between the analytic predictions of the variance of diversity and simulation results. These deviations are likely driven by small deviations in predictions of the second moment of diversity, which are slight for individual community members, but become considerable when terms are summed over hundreds or thousands of community members.

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 (Figure 3a). Similar results were obtained when we repeated our analysis for predicted diversity (Figure 3b). By examining all nine environments we found that despite the dissimilarity in environments, we were able to predict mean richness and diversity in the face of coarse-graining (Figure 3c and d). In contrast, the UNTB failed to predict richness (Figure 3—figure supplement 3). The results of this analysis suggest that the composition of microbial communities remained largely invariant under coarse-graining and that the gamma distribution remained a suitable null model for predicting mean community measures across coarse-grained scales. Identical results were obtained for taxonomic coarse-graining (Figure 3—figure supplement 2).

Figure 3 with 3 supplements see all
The gamma distribution successfully predicted mean richness and diversity under phylogenetic coarse-graining.

(a) The expected richness derived from the gamma distribution (Equation 13a) was capable of predicting richness across phylogenetic coarse-graining scales, as illustrated by data from the human gut. (b) Predictions remained successful across all environments, suggesting that a minimal model of zero interactions was sufficient to predict observed properties of community composition, (c, d) Similarly, predictions of expected diversity (14) also succeeded across coarse-graining scales for all environments. The shade of a color of a given datapoint represents the phylogenetic distance used for coarse-graining, with lighter colors representing finer scales and darker colors representing coarser scales.

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 (Equation 17a). 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 (Figure 4a and b). This lack of predictive success was consistent across environments (Figure 4c and d), implying that a model of independent community members with gamma-distributed abundances was insufficient to capture the variance of measures of biodiversity. A major assumption made in our derivation was that community members are independent, an assumption that is unjustified given that the gamma distribution has been previously shown to be unable to capture the empirical distribution of correlations in the AFDs of community members (Grilli, 2020). 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 (Figure 4a and b) as well as across environments (Figure 4e and 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 (Figure 4—figure supplement 1).

Figure 4 with 1 supplement see all
The gamma distribution only predicts the variance of richness and diversity under phylogenetic coarse-graining when covariance is included.

(a, b) In contrast with the mean, the variance of richness and diversity estimates predicted by the gamma distribution (Equation 17a) failed to capture empirical estimates from the human gut. Predictions are only comparable when empirical estimates of covariance are included in the predictions of the gamma distribution, meaning that dependence among community members is essential to describe the variation in measures of biodiversity across communities. (c, d) This lack of predictive success was constant across environments, (e, f) though the addition of covariance consistently improves our analytic predictions. The color scale used here is identical to the color scale used in Figure 3.

Predicting patterns of richness and diversity between fine and 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 taxonomic and phylogenetic 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 demonstrate the existence of DBD processes 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 Laland et al., 1999; San Roman et al., 2018). 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 (Sg,m) and the richness in the remaining G1 coarse-grained groups: Sg,mαSGg,m, where Gg denotes the exclusion of group g and α is the slope of the relationship. This formulation was proposed by Madi et al., and to ensure commensurability we adopted it here (Madi et al., 2020). Furthermore, keeping with the approach used by Madi et al., fine and coarse-grained measures were compared across increasing taxonomic and phylogenetic scales (e.g. OTU vs. genus, genus vs. family, etc.), (Madi et al., 2020). Using Equation 1, we then defined each of these estimators in terms of the sampling form of the gamma distribution while accounting for sampling

(3a) Sg,m=|g|igP(0|Nm,x¯i,βi)
(3b) SGg,m=(|G|1)gGggP(0|Nm,x¯g,βg)

Similarly, we used Equation 14a to derive predictions for fine and coarse-grained diversity.

(4a) Hg,m=igxln[x]|Nm,x¯g,βg
(4b) HGg,m=gGggxln[x]|Nm,x¯g,βg

By repeating this calculation for all M sites, we obtained vectors of coarse and fine-grained richness estimates for group g from which we inferred the slope of the fine vs. coarse-grained relationship 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. We include a conceptual diagram visualizing this process as a supplement (Figure 5—figure supplement 1).

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 found that they were rarely less than zero (Figure 5a, Figure 5—figure supplement 2a). The few negative slopes inferred from empirical data were extremely small, having absolute values <10−4 and could be treated as zeros. Furthermore, the distribution of slopes follows the same form across environments, suggesting that the slope of the fine vs. coarse-grained relationship 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 (Schluter and Pennell, 2017), is virtually absent in empirical data and cannot be generated from an empirically validated null model of microbial biodiversity.

Figure 5 with 11 supplements see all
The slope of the fine vs. coarse-grained relationship for richness could be predicted by the gamma distribution, but was novel for estimates of diversity.

(a, b) The predictions of the gamma distribution (Equation 3a) successfully reproduced observed fine vs. coarse-grained richness slopes across scales of phylogenetic coarse-graining. (c, d) In contrast, the predictions of the gamma distribution failed to capture diversity slopes (Equation 4a). The color scale used here is identical to the color scale used in Figure 3. Squared Pearson correlation coefficients (ρ2) are computed over all slopes for all taxa across all coarse-graining scales.

However, only observing positive slopes does not necessarily provide support for 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 slope of the fine vs. coarse-grained relationship 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 (Figure 5—figure supplement 3, Figure 5—figure supplement 4). We consolidated these results by taking the mean slope for a given coarse-grained scale, from which we see that the mean slope predicted by the gamma distribution does a reasonable job capturing empirical slopes across environments (Figure 5b, Figure 5—figure supplement 2b). 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 slope of the fine vs. coarse-grained relationship could be sufficiently predicted even though the gamma distribution only succeeded at predicting mean 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.

As a point of comparison, we predicted the slope of the fine vs. coarse-grained relationship for richness using a UNTB model (Madi et al., 2023) (Supporting information). We found that generally, the UNTB slopes deviated from those obtained from empirical data, exhibiting far greater bias and variation around the 1:1 line than what was observed of the SLM (Figure 5—figure supplement 5, Figure 5—figure supplement 6). By examining the mean slope we found that predictions from the UNTB tended to systematically underpredict the observed slope under both taxonomic and phylogenetic coarse-graining (Figure 5—figure supplement 7). Directly comparing the mean relative error of the UNTB predictions to those of the SLM confirms these observations, as the UNTB predictions tended to have larger errors by an order of magnitude (Figure 5—figure supplement 8, Figure 5—figure supplement 9). To summarize, in contrast to the SLM, the UNTB cannot predict the slope of the fine vs. coarse-grained relationship for richness.

While richness is a widespread and versatile estimator that is commonly used in community ecology, neglects considerable 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 was 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 (Figure 5c, Figure 5—figure supplement 2c), 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 (Figure 5d, Figure 5—figure supplement 10, Figure 5—figure supplement 11, Figure 5—figure supplement 2d). 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. Ho et al., 2022), 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 observed fine vs. coarse-grained diversity slopes (Figure 6, Figure 6—figure supplement 1, Figure 6—figure supplement 2, Figure 6—figure supplement 3). 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 (e.g. genus) contribute to measures of biodiversity at the nearest coarse-grained scale (e.g. family), resulting in a positive relationship between measures of diversity at different scales.

Figure 6 with 3 supplements see all
Including correlations allows the gamma distribution to capture observed diversity slopes.

Observed fine vs. coarse-grained diversity slopes could be quantitatively reproduced under phylogenetic coarse-graining by simulating correlated gamma-distributed AFDs at the OTU-level. The color scale used here is identical to the color scale used in Figure 3. Squared Pearson correlation coefficients (ρ2) are computed over all slopes for all taxa across all coarse-graining scales.


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 parameters estimated from data across scales, implying that certain macroecological patterns of microbial communities remained 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 (Whittaker, 1972; Madi et al., 2020). 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. However, introducing correlations between abundance fluctuation distribution permitted the recovery of the slope of the fine vs. coarse-grained diversity relationship.

Our richness results complement past work demonstrating that occupancy, the constituent of richness, is highly dependent on two parameters: sampling depth (i.e. total read count) and the mean abundance of a community member (Grilli, 2020). Our ability to predict the relationship between fine and coarse-grained measures of richness using the gamma distribution, despite our inability to predict the variance of richness, suggest that correlations driving the slope of the fine vs. coarse-grained relationship is primarily driven by the effects of finite sampling. 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 sites (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 the slope of the fine vs. coarse-grained relationship 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 (Madi et al., 2020; Madi et al., 2023). 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 distribution. Contrasting with richness, the predictions of diversity from the gamma distribution were unable to capture fine vs. coarse-grained relationships in empirical data. Given that measures of diversity incorporate information about the richness and evenness of a community (Magurran, 2004), 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 relationships between phylogenetic/taxonomic scales that depend on the evenness of the distribution of abundances.

Macroecological patterns are not imbued with mechanistic explanation (Warren et al., 2022). Rather, the onus is on the investigator to identify plausible mechanisms. Often in ecology this task is made easier by evaluating whether a model lacking a particular mechanism is capable 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 demographic equivalence among community members and community dynamics driven by demographic noise (i.e. the UNTB) (Madi et al., 2020; Alonso and McKane, 2004). 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 reproduces several empirical patterns, making the SLM an appropriate choice for evaluating the novelty of fine vs. coarse-grained relationships (Grilli, 2020; Descheemaeker and de Buyl, 2020). 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 (Williams, 1947; Smith et al., 2014).

In this study, the predictions of the sampling form of the gamma distribution considerably improved when correlations between community members were included. This result suggests that rather than exclusively pointing to niche construction as previously suggested (Madi et al., 2020), 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 (Chesson, 1990; Cui et al., 2021), including quantitatively predicting the distribution of correlation coefficients (Ho et al., 2022), 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 slopes evaluated here have found the existence of positive slopes 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 (Estrela et al., 2022). This mechanism, known as cross-feeding, can be viewed as compatible with the concept of niche construction (San Roman et al., 2018) as well as with the original interpretation of Madi et al., 2020.

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) (Nicoletti et al., 2020; Meshulam et al., 2019). 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 (Sireci et al., 2023).

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 can remain valid for time-series studies where the distribution of microbial abundances remains stationary with respect to time (Faith et al., 2013), as the stationary solution of the SLM has successfully characterized microbial community time-series at both the level of OTUs (Grilli, 2020) and strains (Wolff et al., 2023). Given these past results, we predict that the fine vs. coarse-grained relationship results presented here will remain valid in longitudinal studies where community members fluctuate around a single point with respect to time.

Materials and methods

Data acquisition and processing

Request a detailed protocol

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 (Thompson et al., 2017; Madi et al., 2020). We restricted our analysis to the quality control (QC)-filtered subset of the EMP, which was annotated using the closed-reference database SILVA (Quast et al., 2013) 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 sites from each environment. Summary statistics for each environment are provided (Table 1, Table 2).

Table 1
Summary statistics for the 100 sites randomly selected for each environment.

These statistics reflect the data used for taxonomic coarse-graining, as OTUs lacking taxonomic labels were excluded from taxonomic coarse-graining analyses.

EnvironmentTotal # OTUsMean # OTUsMean # reads
Marine sediment16,1101,393.6740,440.38
Human gut6175599.0932,894.50
Human oral4716537.4544,271.22
Human skin17,9551293.4536,344.13
Freshwater sediment12,2311080.9518,979.56
Microbial mat5087200.248,659.36
Table 2
Summary statistics for the 100 sites randomly selected for each environment.

These statistics reflect the data used for phylogenetic coarse-graining as all OTUs could be used.

EnvironmentTotal # OTUsMean # OTUsMean # Reads
Marine sediment41,3044,167.25106,166.92
Human gut10,190862.7344,031.12
Human oral7062614.9746,104.84
Human skin29,4481817.1248,285.68
Freshwater sediment33,1933,569.5965,582.56
Microbial mat11,869431.4223,216.02

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, as the process of sampling without a replacement would bias the sampling distribution for rare community members (Madi et al., 2020).

Coarse-graining protocol

Request a detailed 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 (Huerta-Cepas et al., 2016).

Deriving biodiversity measure predictions

Request a detailed protocol

While the gamma distribution as the stationary solution of the SLM and the sampling form of the gamma distribution have been previously derived (Grilli, 2020), 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

(5) dxidt=xiτi(1xiKi)Self-limiting growth+στiτixiη(t)Environmental noise

Here τi, Ki, 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 (Gardiner, 2009). The dependence of η(t) at time t on an earlier time η(t) is defined as η(t)η(t)=δ(tt)(Gardiner, 2009). 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 by Madi et al., using an ensemble of sites that belong to the same type of ecosystem rather than the time series of a single site (Madi et al., 2020), 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. George and O’Dwyer, 2022). 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 (Grilli, 2020; Wolff et al., 2023; Descheemaeker and de Buyl, 2020).

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 (Alonso and McKane, 2004; Madi et al., 2020). For the convenience of the reader the predicted richness using the form of the UNTB relevant to this study has been rederived (Supporting information).

The stationary distribution of the SLM can be derived using the Itô ↔ Fokker-Planck equivalence and solving for the stationary solution (Grilli, 2020; Engen and Lande, 1996), resulting in the gamma-distributed AFD. Through the SLM, we can define the mean relative abundance and its squared inverse coefficient of variation as x¯i=Ki(1στi2) and βi=2στiστi, respectively. These are parameters that were estimated from the empirical data and were used below to obtain predictions. Using these definitions and the stationary distribution of Equation 5, we obtained the gamma distribution

(6) P(xi|x¯i,βi)=1Γ(βi)(βix¯i)βiexp[xiβix¯i]xiβi1

When we sequence microbial communities, one obtains read counts rather than actual abundances. Therefore, it is necessary to account for the reality of sampling when we apply 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 (N1) and the typical relative abundance of a community member is much smaller than one (xi1), the binomial can be approximated as a Poisson sampling process with the following probability of sampling n reads

(7) P(n|N,xi)=(Nxi)neNxin!

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 x¯i and βi, the parameters we estimate from the data. This distribution can be obtained by solving the convolution of the Poisson and the gamma distribution (Grilli, 2020). The resulting distribution can be considered a negative binomial distribution if sites have identical sampling depths (Fisher, 1941). Using this distribution, we calculated the probability of obtaining nm reads out of a total sampling depth of Nm for the ith OTU in sample m as

(8a) P(nm|Nm,x¯i,βi)=0P(xi|x¯i,βi)P(nnm|Nm,xi)dxi
(8b) =Γ(βi+nm)nm!Γ(βi)(x¯iNmβi+x¯iNm)nm(βiβi+x¯iNm)βi

This distribution requires two parameters that can be estimated from the data (x¯i and βi) and one parameter that is known (total number of reads, Nm). This equation will be used 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

(9) oi=1MmM(1P(0|Nm,x¯i,βi))

And the second moment of occupancy as

(10) oi2=1MmM(1P(0|Nm,x¯i,βi))2

from which we defined the predicted variance of occupancy

(11) Var(oi)=oi2oi2

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

(12) ε=|Obs.Pred.Obs.|

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

(13a) S=i=1Stotaloi
(13b) =1Mm=1Mi=1Sobs(1P(0|Nm,x¯i,βi))
(13c) =Stotal1Mm=1Mi=1StotalP(0|Nm,x¯i,βi)

where Stotal is the total number of observed community members. Similarly, we derived the expected value of Shannon’s diversity (Magurran, 2004).

(14a) H=1Mm=1MHm
(14b) =1Mm=1Mi=1Sobsxln[x]|Nm,x¯i,βi
(14c) =1Mm=1Mi=1Stotal0NmnNmln[nNm]P(n|Nm,x¯i,βi)dn

In physics parlance, these predictions neglect interactions between community members, also known as mean-field predictions. We then calculated the mean-field prediction of Equation 13a from empirical data. However, there is no known analytic solution for the integral inside the sum of Equation 14a. 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 derived the second moments of richness and diversity.

(15a) S2=(iStotaloi)2
(15b) =1MmM(iStotaloi,m)2
(15c) =1MmMi,joi,moj,mdnP(n|N,x¯i,βi)δN,Nm=1
(15d) =dnP(n|N,x¯i,βi)i,j1MmMδN,Nmoi,moj,m
(15e) =i=1StotaldnP(n|N,x¯i,βi)oi2δN,Nm+ijdnP(n|N,x¯i,βi)oiojδN,Nm
(15f) =i=1Stotaloi2|Nm,x¯i,βi+ijoi|Nm,x¯i,βioj|Nm,x¯j,βj

where δi,j is the Kronecker delta.

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

(16a) H2=(iStotalxiln[xi])2
(16b) =1MmM(iStotalxiln[xi])2
(16c) =1MmMi,j(xiln[xi])(xjln[xj])dnP(n|N,x¯i,βi)δN,Nm
(16d) =1MmMi=1StotaldnP(n|N,x¯i,βi)δN,Nm(xiln[xi])2
(16e) +1MmMijdnP(n|N,x¯i,βi)δN,Nm(xiln[xi])(xjln[xj])
(16f) =1MmMi=1Stotal(xln[x])2|Nm,x¯i,βi
(16g) +1MijmMxln[x]|Nm,x¯i,βixln[x]|Nm,x¯j,βj

Where the expected value of the second moment of the diversity term is defined as (xln[x])2|Nm,x¯s,βs=0Nm(nNmln[nNm])2 P(n|Nm,x¯s,βs)dn. From which we obtained the expected value of the variance

(17a) Var(S)=S2S2
(17b) Var(H)=H2H2

We predicted the mean and variance of richness and diversity separately at each coarse-grained scale. Specifically, we coarse-grained the empirical data, estimate x¯s and βs for each coarse-grained community member, and use these estimates to obtain a prediction for each measure of biodiversity.

It is worth noting why the above functions constitute predictions. To obtain values that we can compare with empirical data we estimated the mean and variance of relative abundance across sites for each community member at a given scale. These parameters were used to obtain the expected value of a community-level measure (e.g. richness) using a function. These functions were derived under the assumption that a given probability distribution (i.e. the gamma) provided an appropriate description of the distribution of relative abundances across sites. We then compared the expected value of a community-level measure to the mean value from empirical data and assessed the similarity between the two values.

Fine vs. coarse-grained relationship slope inference

Request a detailed protocol

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 Equation 3a or Equation 4a within a given coarse-grained group and 3b or Equation 4b among the remaining groups. This ‘leave-one-out’ procedure was originally implemented by 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.) (Madi et al., 2020). 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

Request a detailed protocol

Correlated gamma-distributed AFDs were simulated by performing inverse transform sampling. For each environment with M sites, an M×Sobs matrix Z was generated from the standard Gaussian distribution using the empirical Sobs×Sobs correlation matrix calculated from relative abundances. The cumulative distribution U=Φ(Z)Gaus. was calculated and a matrix of the abundances of community members across sites was obtained using the point percentile function of the gamma distribution and the empirical distribution of mean relative abundances and the squared inverse coefficient of variation of abundances: x¯=x¯1,x¯2,,x¯Sobs, β=β1,β2,,βSobs. To simulate the process of sampling, each community of the resulting M×Sobs matrix of true relative abundances X=Φ(U)Gamma1 was sampled using a multinomial distribution with the empirical distribution of total read counts.

Appendix 1

Supporting information: Investigating macroecological patterns in coarse-grained microbial communities using the stochastic logistic model of growth

UNTB richness predictions

Below we rederive a prediction for richness using the form of the UNTB used by Madi et al., as a point of comparison to the SLM predictions derived in the main manuscript (Madi et al., 2020; Alonso and McKane, 2004). When the size of a metacommunity tends towards an asymptotic limit, the stationary distribution for community members of relative abundance x approaches the following continuous distribution (Vallade and Houchmandzadeh, 2003).

(S1) P(x|θ)dx=θx(1x)θ1dx

where θ is Hubbell’s biodiversity parameter (also known as Fisher’s α Fisher et al., 1943).

Using this distribution, we can obtain an expression for the expected number of community members with n sampled individuals out of a total sample size of N.

(S2) S(n|N,m,θ)=θ01P(n|N,m,x)(1x)θ1xdx

where P(n;N,m,x) is the probability of sampling n individuals of relative abundance x given a total sample size N

(S3) P(n|N,m,x)=(Nn)Γ(n+γx)Γ(γx)Γ(N+γ(1x)n)Γ(γ(1x))Γ(γ)Γ(γ+N)

where γ=m(N1)1m. The function Equation S2 is known as the migration-limited zero-sum multinomial distribution (ZSM) (Alonso and McKane, 2004). As m1, Equation S2 approaches a limiting form known as the metacommunity zero-sum multinomial distribution (mZSM). The process of sampling community members under the mZSM can be represented as a binomial distribution.

(S4) S(n|N,θ)=θ01xN(1x)Nn(1x)θ1xdx

Similar to our analysis using the SLM, the binomial can be approximated as a Poisson distribution.

(S5) S(n|N,θ)=θ01exN(xN)nn!(1x)θ1xdx

The resulting integral can be obtained by using the change of variable y=xN and rearranging terms, then approximating the upper limit of integration as infinity (since N1).

(S6a) S(n|N,θ)=θ01ey(y)nn!(1yN)θ1NydyN
(S6b) =θn01eyyn1(n1)!Gamma distribution(1yN)θ1dy
(S6c) θn0eyyn1(n1)!Gamma distribution(1yN)θ1dy
(S6d) =θn(1YN)θ1

By rearranging terms, we obtain a gamma distribution with shape parameter n and rate parameter 1. Because we integrated over a product with a gamma distribution, the variable Y is a gamma-distributed random variable. We can then expand term in the integral using a Taylor series around Y=n and by noticing that Y=n and Y2=n2+n under a gamma distribution.

(S7) S(n|N,θ)=θn(1nN)θ1+12θ(θ1)(θ2)N2(1nN)θ3+O(N3)

We can then predict the richness of a community by summing the abundances from 1 to N.

(S8) S(N,θ)=n=1NS(n|N,θ)

This quantity represents the total observed richness of a sample from a panmictic infinite metacommunity after accounting for sampling. Predictions of mean richness over M sites can then be calculated as

(S9) S(θ)=1Mm=1MS(Nm,θ)
Simulating fine vs. coarse-grained richness slopes under the UNTB

We followed the procedure in Madi et al., to obtain fine vs. coarse-grained slopes for richness so that they could be compared to predictions obtained from the SLM (Madi et al., 2020). We simulated SADs according to the mZSM model outlined above using the rmzsm() function from the R package sads v0.4.2. We simulated 100 SADs using the empirical distribution of total read counts and the total number of observed OTUs. We set the biodiversity parameter θ=50 for all environments. The SADs returned by rmzsm() contain no zeros, meaning that values of richness are identical for all UNTB SADs. In order to introduce zeros so that richness estimates could vary, we followed the procedure used in Madi et al., where each simulated SAD was rarefied to 5000 individuals. We repeated this rarefaction procedure on the empirical SADs. We then performed taxonomic and phylogenetic coarse-graining and fine vs. coarse-grained slope inference using the procedure described in the Materials and methods.

Data availability

All sequencing data used in this study was obtained from the Earth Microbiome Project (URL: Processed data used to perform the analyses in this study are available on Zenodo, DOI: All code written for this study is available on GitHub under a GNU General Public License: (copy archived at Shoemaker, 2023b).

The following data sets were generated
    1. Shoemaker WR
    (2023) Zenodo
    Macroecological patterns in coarse-grained microbial communities.
The following previously published data sets were used


  1. Book
    1. Gardiner CW
    Stochastic Methods: A Handbook for the Natural and Social Sciences
    Berlin Heidelberg: Springer.
  2. Book
    1. Magurran AE
    Measuring Biological Diversity
    Malden, Ma: Blackwell Pub.

Peer review

Reviewer #1 (Public Review):

Shoemaker and Grilli analyze publicly available sequencing data to quantify how the microbial diversity of ecosystems changes with the taxonomic scale considered (e.g., diversity of genera vs diversity of families). This study builds directly on Grilli's 2020 paper which used this data to show that for many different microbial species, the distribution of abundances of the species across sampling sites belongs to a simple one-parameter family of gamma distributions. In this work, they show that the gamma distribution also describes the distribution of abundances of higher taxonomic levels. The distribution now requires two parameters, but the second parameter can be approximately derived by treating the distributions of lower-level taxonomic units as being independent. The difference between the species-level result and the result at higher taxonomic levels suggests that in some sense microbial species are ecologically meaningful units.

While the higher-level taxon abundance distributions can be well-approximated assuming independence of the constituent species, this approach substantially underestimates variation in community richness and diversity among sampling sites. Much of this extra variability appears to be driven by variability in sample size across sites. It is not clear to me how much this variation in sample size is itself due to variation in sampling effort versus variation in overall microbial densities. This variation in sample size also produces correlations between taxon richness at lower and higher taxonomic levels. For instance, sites with large samples are likely to have both many species within a genus and many genera. The authors also consider taxon diversity (Shannon index, i.e. entropy), which is constructed from frequencies and is therefore less sensitive to sample size. In this case, correlations between diversity across taxonomic scales instead appear to depend on the idiosyncratic correlations among species abundances.

This paper's results are presented in a fairly terse manner, even when they are describing summary statistics that require a lot of thought to interpret. I don't think it would make sense to try to understand it without having first worked through the 2020 paper. But everyone interested in a general understanding of microbial ecology should read the 2020 paper, and once one has done that, this paper is worth reading as well simply for seeing how the major pattern in that paper shifts as one moves up in taxonomic scale.

Reviewer #3 (Public Review):


In this research advance, the authors purport to show that the unified neutral theory of biodiversity (UNTB) is not a suitable null model for exploring the relationship between macroecological quantities, and additionally that the stochastic logistic growth model (SLM) is a viable replacement. They do this by citing other studies where UNTB was unable to capture individual macroecological quantities, and then demonstrating SLM's strength at predicting the same diversity metrics. They extend this analysis to show SLM's modeling capability at multiple scales of coarse graining, in addition to its failures at predicting these metrics' variances. Finally, authors conduct a similar analysis to Madi et al. (2020) by investigating the relationship between diversity measures within a group and across coarse-grained groups (e.g. genera diversity in one family compared to diversity of families). The authors show that choosing SLM as a null model reveals some previously reported relationships to be no longer "novel", in the sense that the patterns can be adequately captured by the null model. Authors also show that relationships not captured by the null model can be recovered by adding correlations, suggesting interactions are the driving force behind them.


1. Authors make a strong argument that UNTB is not a good null model of macroecological observables and especially relationships between them. Authors convincingly argue that a SLM is a better null since the gamma distribution it predicts is a better description of the empirical Abundance Fluctuation Distributions (AFD).

2. Authors show that the gamma distribution predicted by SLM is a good fit for the AFD's at many different scales of coarse graining, not just the OTU level as was previously demonstrated. Authors show the same distribution predicted the mean diversity and richness at all scales of coarse graining.

3. Authors convincingly demonstrate how SLM can be used to test the relevance of interactions to macroecological relationships.


This reviewer's concerns were convincingly addressed by the revisions.

Overall Impact

The authors present a convincing argument for the use of SLM as a better non-interacting null model for macroecological quantities and relationships.

Author response

The following is the authors’ response to the original reviews.

Reviewer #1:

In no particular order:

1. In Figs S3 and S4, can they also show gamma fit? (or rather corrected fit accounting for abundance conditioning?) The shapes look different, especially for the microbial mat.

Author response: We have added gamma distribution fits to the rescaled AFD plots (Figs. S3, S4).

1. Lines 170-176 seem like they should come before lines 164-166.

Author response: In lines 166-170 we discuss empirical patterns in the data that motivate the introduction of the SLM as a model in lines 170-175. We have clarified these points in the revision.

1. The wiggles in the gamma predictions in the occupancy-abundance plots are because occupancy depends not only on abundance but also on the shape parameter, right? Probably good to write a sentence or two explaining what's going on here.

Author response: We agree with the reviewer that the variation in the prediction could be in-part driven by variation in the shape parameter across community members. We now include this observation in our revision (lines 209-211).

1. In the predicted vs observed occupancy plots, it would be nice to add curves showing predicted standard deviation or similar to give a sense of how well the model is predicting the variability.

Author response: In the revised manuscript we now include predictions for the variance of occupancy using the gamma distribution under both taxonomic and phylogenetic coarse-graining (Fig. S9; S10; lines 211-214).

1. Covariance between sister groups: Figs S9 and S10 look very nice, but it's hard to see much because they're log-log plots over multiple decades, while even a several-fold difference from y = x would indicate a strong effect of correlations. It would be clearer if the y-axis showed the ratio of the coarsegrained variance to the sum of OTU variances and we were looking at how well it fit y = 1.

Author response: We have included these plots in the revision (Fig. S14, S15).

1. If the sum of gammas can be well-approximated by a gamma, does that mean that the gamma is just a fairly flexible distribution and we shouldn't take the quality of the gamma fits in general as a very specific indication of what's going on?

Author response: While the sum of random variables that are drawn from gamma distributions with different parameters is often well-approximated by another gamma, this does not tell us why the gamma distribution holds for microbial communities at the finest-grain level (i.e.,OTUs/ASVs). At present, the best explanation is that the gamma is a stationary distribution for certain stochastic differential equations which have ecological interpretations (Grilli, 2020; Shoemaker et al., 2023). Furthermore, alternative two-parameter distributions have been tested alongside the gamma and have done a comparatively poor job capturing observed macroecological patterns (Grilli, 2020). These results suggest that the utility of the gamma distribution is not simply an outcome of its flexible nature, it succeeds because it has captured core ecological properties of microbial communities. In the case of the SLM, gamma-like distributions arise when a community member is subject to self-limiting growth and environmental noise. On the other hand, the stability of the gamma distribution might explain why it can be detected as shape of the AFD, as it does not fade out across coarse-graining level.

1. What's going on with the variance of diversity in Fig S12? Does this suggest that some of the problem in Figure 4 could be with the analytic approximation rather than the model? I had a hard time understanding the part of the Methods explaining the simulation details (lines 587-597). It would be worth expanding this. Is there some way to explain how the correlations were simulated in terms of the SLM, e.g., correlations in the noise term across OTUs?

Author response: We believe that deviations in the variance of diversity in Fig. S16g,h are driven by small deviations in our predictions of the second moment (xln(x)|Nm,x¯i,βi2) (Eq. S16). Alone these predictions are slight, but their effects become noticeable when summed over hundreds or thousands of taxa. We have included this observation in the revised manuscript (lines 268-271). However, this deviation pales in comparison with the magnitude of covariance in the empirical data, suggesting that our inability to predict the variance of richness and diversity is primarily driven by our assumption of statistical independence.

Regarding the source of the correlations, under the SLM correlations in abundances can be introduced either by adding deterministic interaction terms or through correlated environmental noise. Determining which of these two options drives empirical correlations is an active area of research (e.g., Camacho-Mateu et al., 2023). For the purpose of this study, we remain agnostic on the cause of the correlations, optioning to instead emphasize that that the inclusion of correlations is necessary to reproduce observed slopes of the fine vs. coarse-grained relationship for diversity.

1. In Figure 5ab, is the idea that the correlation in richness is primarily driven by the number of samples from the environment? Line 390 seems to say so, but it would be good to make this explicit and put it right in that section of the Results.

Author response: Our results suggest that sampling effort (# reads) plays a larger role in determining the correlations between fine and coarse-grained measures of richness. We now clarify this point in the revised manuscript (lines 429-435).

1. I don't totally understand the contrast in lines 369-372. If fine-scale diversity within one group begets coarse-grained diversity in another group, couldn't that show up as correlations in the AFDs? Or is the argument that only including within-group correlations in AFDs is enough to reproduce the pattern? I'm not sure I see how that could be.

Author response: The term “begets” implies both causation and direction. If we see a positive relationship between diversity estimates at two different scales of observation the causal mechanism cannot be determined solely from correlations between samples obtained once from different sites. So, mechanisms consistent with niche construction/"DBD" can produce correlations, though the existence of correlations do not necessarily imply DBD.

1. The discussion of niche construction on 429-431 doesn't match very well with 440-441. Basically, niche construction is a very broad concept, not a specific one, right?

Author response: In lines 472-576 (formerly 429-431) we discuss how the existence of correlations between fine and coarse-grained scales does not point to a single ecological mechanism. Alternatively stated, observing a non-zero slope does not mean that niche construction is driving the relationship.

In lines 476-487 (formerly 440-441) we discuss how the mechanism of cross-feeding has been shown to generate a positive relationship between fine and coarse-grained measures of diversity. This mechanism can be interpreted as a form of “niche construction”, so it is an instance of a tested ecological mechanism that aligns with the interpretation given in Madi et al. (2020).

1. Isn't (8) just the negative binomial distribution?

Author response: The convolution of the stationary solution of the SLM (i.e., a gamma distribution) and the Poisson limit of a multinomial sampling distribution returns a negative binomial distribution of read counts across hosts if samples have identical sampling depths. We now include this detail in the revision (line 593-595). Note however that if different samples have different sampling depths, the distribution of reads across samples is not a negative binomial.

1. Missing 1/M in (9).

Author response: We have fixed this omission in the revision.

1. Schematic figures illustrating what the different statistics are intuitively capturing would really help this work be understandable to a broader audience, but they'd also be a ton of work.

Author response: Richness and diversity are used in ecology to such an extent that we do not see the benefit of a conceptual diagram. Furthermore, we have included a conceptual diagram about our pipeline in our revision at the request of Reviewer 2 (Fig. S20).

Reviewer #2:

Major Recommendations

If I were reviewing this manuscript for a regular journal, I believe the following issues would be important to address prior to publication.

1. From my reading, the main points of this advance are that

a. SLM models AFDs well at all levels of coarse-graining.

b. This makes SLM a better null-model than UNTB for macroecological relationships.

c. Using SLM on the EMP data, the richness slopes are well explained by SLM but not the diversity slopes. Therefore, any theory that hopes to explain the diversity slopes must include interactions. Argument B appears to be one of the key points yet is missing from the abstract, and should be made clearer. If these aren't the main points the authors intended, then other main points need to be highlighted more.

Author response: In the revision we now explicitly mention argument b in the Abstract.

1. The title should be more specific, so as to better reflect the content. (E.g. "UNTB is not a good null model for macroecological patterns" would seem more appropriate.)

Author response: We would prefer to focus on the success of the SLM rather than the limitations of the UNTB in the title of this work. Therefore, we have modified our title as follows: “Investigating macroecological patterns in coarse-grained microbial communities using the stochastic logistic model of growth”.

1. The manuscript would benefit from a clearer description of exactly what information the SLM retains about the data (perhaps even a cartoon panel in one of the figures). In particular, it is important to be explicit about the number of model parameters.

Author response: The number of model parameters for the gamma AFD are now explicitly stated in the revision (Lines 579-580).

1. The main point of Figures 2-4 seems to be that SLM is good at describing the data (and when it fails it is due to interactions) while UNTB fails to reproduce this behavior, in support of Argument B. This is not clear from the figure descriptions or titles, which focus on SLM's "predictive" power.

Author response: Fig. 2a demonstrates that the gamma distribution predicted by the SLM explains the empirical distribution of abundances. This result provides motivation to predict the fraction of sites harboring a given community member (i.e., occupancy, Fig. 2c) as well as general measures of community composition including mean richness (Fig. 3a,c) and mean diversity (Fig. 3b,d) using parameters estimated from the data (not free parameters).

This success led us to consider whether the gamma distribution could predict the variance of richness and diversity, which it could not because it does not capture covariance between community members (Fig. 4).

In the revision we have identified opportunities to make these points clear throughout the Results.Furthermore, we have added additional detail to the legends of Figs. 2-4.

1. The manuscript would benefit from clarifying the use of "prediction" related to the SLM. Since the gamma distributions predicted by SLM were fit to empirical data, it seems like the agreement between analytic means and empirical means (Fig. 3) is a statement on gamma distributions being a good fit for the AFD's more than SLM predicting richness and diversity. For example, from my reading, it seems like this analysis could be done numerically by shuffling species abundances across environments and seeing whether this changed the mean richness/diversity. I would not call this shuffling test a prediction, since it is more a statement on the relevance of interactions. SLM predicts gamma-distributed AFD's, but those distributions recovering the data they were trained on doesn't seem like a prediction.

Author response: In this manuscript we identified the gamma distribution as an appropriate probability distribution to describe the distribution of relative abundances across samples over a range of coarse-grained scales. Motivated by this result, we performed a separate analysis where at each scale we estimated the mean and variance of relative abundance across sites for each community member. We then used these parameters to obtain the expected value of acommunity-level measure using an equation we derived by assuming that the gamma distribution was appropriate (e.g., richness, Eq. 13). We then compared the expected value of richness to the mean value from empirical data and assessed the similarity between the two values.

The outcome of this procedure constitutes a prediction. While the mean and variance are parameters, estimating them from the empirical data has no connection with the operation of training a distribution on empirical data. We could have derived predictions such as Eq. 13 using any other probability distribution that can be parameterized using the mean and variance (e.g., Gaussian). Such a prediction would likely do a poor job even though it used the same means and variances used for our gamma predictions. This is because the choice of distribution would not have been a good descriptor of the distribution of abundances across hosts.

To better explain this last -- perhaps the most significant -- issue, I'd like to ask the authors if the following recasting would be an accurate reflection of their conclusions, or if something is missing.

1. "Focusing on the empirical relationship observed between diversity slopes by Madi 2020, we ask the question: does explaining these relationships require accounting for species-species correlations? Or could it be reproduced in a noninteracting model?"To address this question, one can perform a randomization test, shuffling abundances to preserve all single-OTU statistics but breaking any correlations. My reading of the authors' results is that (new result 1) the richness relationships would be preserved, while diversity relationships would not be preserved. [Note that this result 1 need not mention either SLM or UNTB.]

Author response: The question of whether correlations between species are necessary to explain the observed slope of the fine vs. coarse-grained relationship was only one component of our research goals. Our first question was whether the SLM would prove to be a more appropriate null for evaluating the novelty of observed slopes. We believe that our results support the conclusion that the SLM is an appropriate null for this question, as it was able to capture observed slopes of the fine vs. coarse-grained relationship for estimates of richness, determining that correlations and the interactions that are ultimately responsible are not necessary to explain this result.

We then find that the SLM as a null model fails to capture observed slopes of the fine vs. coarsegrained relationship for estimates of diversity and simulate the SLM with correlations to return reasonable estimates of the slope. However, here the question about correlations is a direct follow-up from our question about a null model that excludes interactions, so it is unclear how a randomization test would relate to this result.

1. Instead of doing a randomization test (resampling the empirical distribution), one might insist on instead fitting a model to the AFD distributions, and sampling from that distribution rather than the empirical one.

a. If doing it this way, one should of course ensure that the distribution being fit is a good description of the data.

b. UNTB is a bad fit. SLM is a better fit, and in fact (new result 2) continues to be a good empirical fit even at coarse-grained levels.

c. Can make statements on using SLM as a null model for these types of cross-scale relationships. Could try arguing that fitting an SLM model per-OTU (instead of resampling the empirical distribution) could offer some advantage if certain properties could be computed analytically from the fit parameters, instead of averaging over multiple computational rounds of resampling.

Do these two points accurately summarize the manuscript? If so, this presentation avoids the confusion with "prediction". If my summary is missing some important point, the presentation should be revised to clarify the points I appear to have missed.

Author response: In our manuscript we derive predictions from the gamma distribution, the stationary distribution of the SLM, that require parameters estimated from the data (i.e., mean and variance of relative abundance). These parameters are estimated from the data using normal procedures and then plugged into our predictions that assume the appropriateness of the gamma, returning values that are then compared to estimates from empirical data. Our estimation of the mean and variance does not assume that the empirical distribution following a gamma distribution, but the value returned by our function derived from the gamma distribution (e.g., Eq. 13) does make that assumption.

To address the reviewer’s broader comment, we believe that following points summarize our manuscript:

1. The gamma distribution as a stationary solution of the SLM captures macroecological patterns and predicts typical community-level properties (i.e., mean richness and diversity) across phylogenetic and taxonomic scales.

2. The gamma distribution fails to predict variation in community-level properties (i.e., variance of richness and diversity) across phylogenetic and taxonomic scales. This occurs because the SLM is a mean-field model that does not explicitly include interactions between community members.

3. Despite the inability to capture interactions, the gamma distribution succeeds at predicting the fine vs. coarse-grain slope for richness, a pattern that had previously been attributed to community member interactions. This result demonstrates that the novelty of a macroecological pattern hinges on one’s choice of null model.

4. However, the gamma cannot capture the same relationship for diversity. Simulations of the gamma distribution that incorporate correlations between community members are capable of generating reasonable estimates of the slope.

To address the reviewer’s comments regarding the appropriateness fitted gamma distributions, in our revision we have added fitted gamma distributions to plots of AFDs so that the reader can visually assess the ability of the gamma to describe empirical patterns (Fig. S3, S4).

We have also obtained predictions for the slope of the fine vs. coarse-grained relationship for community richness using the same form of UNTB used by Madi et al (2020). In our revised manuscript we establish a procedure to infer the single parameter of this model, generate predictions of richness at fine and coarse-grained scales, and then evaluate whether the UNTB is capable of predicting the slope of the fine vs. coarse-grained relationship for richness (Supplementary Information; Figs. S18, 24-28; lines 277-278; 370-380).

Other/minor comments

1. The manuscript would be improved with more consistent terminology ("fine vs. coarse-grained relationship"/"the relationship" vs. "diversity slope"). Also, many readers may be used to OTUs referring to the rather fine level of description, as opposed to any chosen level; and could interpret indexing over groups as being in contrast with indexing over OTU's (coarse vs fine). The authors' use is perfectly correct, but keeping a consistent terminology would help.

Author response: We have revised our manuscript to specify the “slope” as the “slope of the fine vs. coarse-grained relationship” (e.g., Line 318). We also specify in the Results and in the Methods that we use “fine” and “coarse” as relative terms, keeping with the sliding-scale approach used in Madi et al (2020).

1. While I appreciate this "slope" is something borrowed from other work, the clarity of the paper might benefit from a cartoon of how one goes from the raw data to the slopes at a particular coarse-graining level. (Optional).

Author response: We had added a conceptual diagram to the revision (Fig. S20).

1. The text often colloquially references "the gamma," "predictions of the gamma," etc. This phrasing comes across as sloppy, and the manuscript would be improved by being more specific.

Author response: We now specify “gamma” as the “gamma distribution” throughout the manuscript.

1. Equation 6 appears to be missing some subscripts on the x terms (included on the left of the equation).

Author response: We thank the reviewer for noticing this error and we have corrected it in the revision.

1. In "Simulating communities of correlated...AFDs", the acronym SAD is not defined.

Author response: We thank the reviewer for noticing this error and we have corrected it in the revision.

1. In Figure 2:

a. Invariant is probably the wrong word for the title, since all the AFD's were rescaled by mean and variance before being compared. Data does support that the gamma distributions are good at describing the AFD's, but as stated in the description it's the general shape that is preserved, not the distribution itself.

Author response: When we mention the invariance of the AFD we now specify that we mean that the shape of the distribution remained qualitatively invariant.

b. I'd recommend changing the color coding to something with more contrast, since currently it's impossible to assess the claim that the shape of the distribution collapses.

Author response: Our coarse-graining procedure is a sequential operation that has no intuitive point that would suggest the use of a contrasting colormap (e.g., if our scale ranged from -1 to 1 then there would be a natural point of contrast at zero).

c. The legend is missing relevant technical details: How many OTU's were used to make plot a? How many samples?

Author response: The number of samples was listed in the Materials and Methods (line 523). In the revision we now include a table with the average and total number of OTUs as well as the average number of reads for each environment (Table S1, S2).

d. In plot b, is the mean relative abundance referring to "mean abundance when observed" or "mean across all samples"?

Author response: The mean relative abundance is the mean abundance across all sites (line 204) and in the legend of Fig. 2.

e. Since one argument here is that SLM fits these distributions better than UNTB, if possible it would be nice to see UNTB's failed fits here.

Author response: A major feature of the UNTB is that the demographic parameters of community members are indistinguishable. Under the SLM, the variation in the mean relative abundance we observe suggests that the carrying capacities of community members vary over multiple orders of magnitude, a result that is incompatible with most forms of the UNTB (x-axis of Fig. 2b). We now mention this point in the revised manuscript (lines 110; 229; 455-471).

1. In Figure 3:

a. It is not clear how coarse-graining is included in model fitting. The "Deriving biodiversity measure predictions" section would benefit from including how coarse-graining is incorporated.

Author response: We predict measures of biodiversity separately at each coarse-grained scale. We now clarify this detail in the revised manuscript (Lines 624-627).

b. Reference Shannon Diversity in Methods.

Author response: We now cite Shannon’s diversity.

c. What is the blue/white color coding in plots a & c? It doesn't have any color key.

Author response: Figs. 3-6 use a uniform light-to-dark scale for all environments, with each environment having its own color. For example, Fig. 3a contains data from the human gut microbiome. Human gut data were assigned the color aquamarine, so the shade of aquamarine for a given datapoint in Fig. 3a indicates the phylogenetic scale.

In the revision we now clarify the colorscale in the legend of Fig. 3 and specify that the same scale is used in all subsequent figure legends.

d. Re: earlier comments, why is richness considered a prediction? (Am I correct in my interpretation that panel b is almost a tautology - counting the number of zeros in the matrix either by rows or by columns - whereas panel d is nontrivial?)

Author response: Mean richness as a measure of biodiversity depends on the fraction of sites where a given community member is present (i.e., occupancy). The mean relative abundance of a community member and its variation across sites (beta) is clearly related to occupancy, but those two statistics do not give you a prediction of occupancy. Obtaining a prediction of occupancy and, subsequently, richness, requires 1) a probability distribution of abundances (i.e., the gamma) and 2) a probability distribution of sampling (i.e., the Poisson). Using these two pieces of information, we derived a prediction for mean richness (Eq. 13). We then compare the value of richness obtained by plugging in the mean relative abundances, betas, and known number of reads to the observed mean richness obtained from the data.

e. The lettering of subplots in Figure 3 is not consistent with Figure 4. Figure 3 subplots are also cited incorrectly in paragraph two on page six (lines 251-254).

Author response: We thank the reviewer for noticing the error and we have corrected it in the revision.

f. Again, if possible show UNTB predictions in plots a & c.

Author response: In our revised manuscript we provide extensive descriptions and predictions of mean richness and the slope of the fine vs. coarse-grained relationship for richness using the form of the UNTB used in Madi et al. (2020; Figs. S18, S24 - S29; lines 277-282; 370-380). We then compare the error of these slope predictions to those obtained from the SLM, finding that the SLM generally outperforms UNTB (Figs. S27-S29).

1. In Figure 4:

a. What are the color codings in plots a & b?

Author response: The color scale used in Fig. 4 is identical to the color scale used in Fig. 3. This detail is now specified in the legend of Fig. 4.

b. What are the two lines of empirical data in plots a & b, and why is one of them dashed?

Author response: We now specify what the two lines mean in the key within the figure.

c. Same comment as earlier on predictions and richness.

Author response: We now specify what the two lines mean in the key within the figure.

1. In Figure 5:

a. It wasn't clear to me in the manuscript how the authors generated these plots from the raw data. The manuscript would benefit from a clear cartoon/description of the data pipeline, from raw data to empirical (and analytic) slopes.

Author response: We have added a conceptual diagram to the revised manuscript (Fig. S20).

b. Make the figure title more descriptive to better connect it to the figure's objective (the richness slopes relationship is not novel, but the diversity slopes relationship is).

Author response: We have revised the figure title.


Camacho-Mateu, J., Lampo, A., Sireci, M., Muñoz, M. Á., & Cuesta, J. A. (2023). Species interactions reproduce abundance correlations patterns in microbial communities(arXiv:2305.19154). arXiv.

Grilli, J. (2020). Macroecological laws describe variation and diversity in microbial communities. Nature Communications, 11(1), 4743.

Madi, N., Vos, M., Murall, C. L., Legendre, P., & Shapiro, B. J. (2020). Does diversity beget diversity in microbiomes? eLife, 9, e58999.

Shoemaker, W. R., Sánchez, Á., & Grilli, J. (2023). Macroecological laws in experimental microbial systems (p. 2023.07.24.550281). bioRxiv.

Article and author information

Author details

  1. William R Shoemaker

    Quantitative Life Sciences, The Abdus Salam International Centre for Theoretical Physics (ICTP), Trieste, Italy
    Conceptualization, Resources, Data curation, Software, Formal analysis, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0111-4838
  2. Jacopo Grilli

    Quantitative Life Sciences, The Abdus Salam International Centre for Theoretical Physics (ICTP), Trieste, Italy
    Conceptualization, Formal analysis, Supervision, Methodology, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared


National Science Foundation (2010885)

  • William R Shoemaker

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.


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

Senior Editor

  1. Meredith C Schuman, University of Zurich, Switzerland

Reviewing Editor

  1. Bernhard Schmid, University of Zurich, Switzerland

Version history

  1. Preprint posted: May 24, 2023 (view preprint)
  2. Sent for peer review: June 19, 2023
  3. Preprint posted: August 21, 2023 (view preprint)
  4. Preprint posted: December 21, 2023 (view preprint)
  5. Version of Record published: January 22, 2024 (version 1)

Cite all versions

You can cite all versions using the DOI This DOI represents all versions, and will always resolve to the latest one.


© 2023, Shoemaker and Grilli

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.


  • 234
    Page views
  • 27
  • 0

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. William R Shoemaker
  2. Jacopo Grilli
Investigating macroecological patterns in coarse-grained microbial communities using the stochastic logistic model of growth
eLife 12:RP89650.

Share this article

Further reading

    1. Ecology
    Congnan Sun, Yoel Hassin ... Yossi Yovel
    Research Article

    Covid-19 lockdowns provided ecologists with a rare opportunity to examine how animals behave when humans are absent. Indeed many studies reported various effects of lockdowns on animal activity, especially in urban areas and other human-dominated habitats. We explored how Covid-19 lockdowns in Israel have influenced bird activity in an urban environment by using continuous acoustic recordings to monitor three common bird species that differ in their level of adaptation to the urban ecosystem: (1) the hooded crow, an urban exploiter, which depends heavily on anthropogenic resources; (2) the rose-ringed parakeet, an invasive alien species that has adapted to exploit human resources; and (3) the graceful prinia, an urban adapter, which is relatively shy of humans and can be found in urban habitats with shrubs and prairies. Acoustic recordings provided continuous monitoring of bird activity without an effect of the observer on the animal. We performed dense sampling of a 1.3 square km area in northern Tel-Aviv by placing 17 recorders for more than a month in different micro-habitats within this region including roads, residential areas and urban parks. We monitored both lockdown and no-lockdown periods. We portray a complex dynamic system where the activity of specific bird species depended on many environmental parameters and decreases or increases in a habitat-dependent manner during lockdown. Specifically, urban exploiter species decreased their activity in most urban habitats during lockdown, while human adapter species increased their activity during lockdown especially in parks where humans were absent. Our results also demonstrate the value of different habitats within urban environments for animal activity, specifically highlighting the importance of urban parks. These species- and habitat-specific changes in activity might explain the contradicting results reported by others who have not performed a habitat specific analysis.

    1. Ecology
    2. Plant Biology
    Daniel Fuks, Yoel Melamed ... Ehud Weiss
    Research Article Updated

    Global agro-biodiversity has resulted from processes of plant migration and agricultural adoption. Although critically affecting current diversity, crop diffusion from Classical antiquity to the Middle Ages is poorly researched, overshadowed by studies on that of prehistoric periods. A new archaeobotanical dataset from three Negev Highland desert sites demonstrates the first millennium CE’s significance for long-term agricultural change in Southwest Asia. This enables evaluation of the ‘Islamic Green Revolution (IGR)’ thesis compared to ‘Roman Agricultural Diffusion (RAD)’, and both versus crop diffusion during and since the Neolithic. Among the findings, some of the earliest aubergine (Solanum melongena) seeds in the Levant represent the proposed IGR. Several other identified economic plants, including two unprecedented in Levantine archaeobotany—jujube (Ziziphus jujuba/mauritiana) and white lupine (Lupinus albus)—implicate RAD as the greater force for crop migrations. Altogether the evidence supports a gradualist model for Holocene-wide crop diffusion, within which the first millennium CE contributed more to global agricultural diversity than any earlier period.