Robust estimation of microbial diversity in theory and in practice

Abstract

Quantifying diversity is of central importance for the study of structure, function and evolution of microbial communities. The estimation of microbial diversity has received renewed attention with the advent of large-scale metagenomic studies. Here, we consider what the diversity observed in a sample tells us about the diversity of the community being sampled. First, we argue that one cannot reliably estimate the absolute and relative number of microbial species present in a community without making unsupported assumptions about species abundance distributions. The reason for this is that sample data do not contain information about the number of rare species in the tail of species abundance distributions. We illustrate the difficulty in comparing species richness estimates by applying Chao's estimator of species richness to a set of in silico communities: they are ranked incorrectly in the presence of large numbers of rare species. Next, we extend our analysis to a general family of diversity metrics ('Hill diversities'), and construct lower and upper estimates of diversity values consistent with the sample data. The theory generalizes Chao's estimator, which we retrieve as the lower estimate of species richness. We show that Shannon and Simpson diversity can be robustly estimated for the in silico communities. We analyze nine metagenomic data sets from a wide range of environments, and show that our findings are relevant for empirically-sampled communities. Hence, we recommend the use of Shannon and Simpson diversity rather than species richness in efforts to quantify and compare microbial diversity.

Figure 1.

Empirical sample data are consistent with very different communities. We consider the abundance data of a sample taken from a bacterial soil community (sample ‘Brazil' in (Roesch et al., 2007)). The sample consists of 26079 individuals belonging to 2880 species. We tried to reconstruct the community from which the sample was taken. Panels a–c show the rank-abundance curve of three such reconstructed communities. The first community (panel a, in red) has 10⁴ species; the second community (panel b, in blue) has 10⁵ species; the third community (panel c, in green) has 10⁶ species. For each of the three reconstructions the community rank-abundance curve is an extension of the sample rank-abundance curve (in black). We claim that each of the three reconstructed communities is compatible with the sample data. This can be seen from the rarefaction curves in panel d: the rarefaction curve for the sample data (black line) coincides with the rarefaction curves for the reconstructed communities (red line with squares for community in panel a, blue line with x-marks for community in panel b, and green line with diamonds for community in panel c). Because the sample data are consistent with very different values of the community richness, the community richness cannot be estimated from the sample data. The colour reproduction of this figure is available on the ISME Journal online.

Figure 2

Estimated species richness does not rank correctly communities. We generated three community abundance distributions, the rank-abundance curves of which are shown in panel a. Community C1 (red) has the smallest number of species; community C3 (green) has the largest number of species. The rarefaction curves of the three communities up to sample size 2 10⁴ are shown in panel b. Based on the rarefaction data, one would conclude that community C1 is the most diverse and community C3 the least diverse. Hence, the ranking of the communities according to their observed diversity is inverted compared to the ranking according to their true diversity. This observation is confirmed when applying Chao's estimator to sample data. Community C1 is estimated to have 10 times more species than community C3, whereas in reality community C1 has 20 times less species than community C3. See Supplementary Table S1 for the numerical data of the communities. The colour reproduction of this figure is available on the ISME Journal online.

Figure 3

Extrapolating the rarefaction curve. The Hill diversity estimators and are based on reconstructions of the rarefaction curve S_m from sample data. For a sample of size M, the rarefaction curve S_m for m⩽M can be estimated by subsampling (red full line). If the sample size M is large, the estimator has small uncertainty. The rarefaction curve S_m for m>M can be estimated by extrapolating the sample data beyond the sample size M. Different extrapolation scenarios are compatible with the sample data. We consider two extreme scenarios (dashed lines). A lower estimate is obtained by assuming that unobserved species are approximately as rare as the rarest observed species. An upper estimate is obtained by assuming that unobserved species are represented in the community by one individual. The difference between the two extremes quantifies the uncertainty of the extrapolation, shown as the shaded region. The uncertainty increases rapidly for m>>M.

Figure 4

Estimated Hill diversities for in silico communities. We generated samples from a community with power-law abundance distribution (S=10⁶, z=2) and evaluated the estimators and for the Hill diversity D_α. We consider three sample sizes M (in columns: M=10², 10⁴, 10⁶) and three community sizes N (in rows: N=10¹⁰, 10¹⁵, 10²⁰). The shaded range between and indicates the estimation uncertainty. The true Hill diversity D_α of the community is plotted in black. The Hill diversities between α=1 (Shannon) and α=2 (Simpson) are correctly estimated even for small sample size M. The estimates of Hill diversities less than α=1, including α=0 (species richness), are characterized by large uncertainty.

Figure 5

Estimated Hill diversities for natural microbial communities. We observe the same behavior as for the in silico generated data sets of Figure 4: for α⩾1 the Hill diversity D_α can be estimated accurately; for α<1 the estimation of the Hill diversity D_α has large uncertainty. We used the same data sets as Quince et al., (2008): a seawater bacterial sample from the upper ocean (Rusch et al., 2007), soil bacterial samples at four locations: Brazil, Florida, Illinois and Canada (Roesch et al., 2007), and seawater samples from deep-sea vents at two locations: FS312 and FS396, separated into bacteria and archaea (Huber et al., 2007). The community size was set to N=10¹⁵ for illustration; results are robust to changes in community size (see Supplementary Figure S4).