Alpha-, beta-, and gamma-diversity of bacteria varies across habitats

1Bacteria are the most diverse organisms on the planet. Bacterial richness and composition influences ecosystem functioning, whether in host-associated communities, soils, or oceans. Nevertheless, we have yet to answer a number of basic questions about bacterial diversity, including “Which habitats contain the highest diversity of bacteria?” More broadly, evaluating geographic patterns in biodiversity across habitats and spatial scales can illuminate the processes influencing and consequences of biodiversity.

2While many studies document spatial patterns of bacterial diversity, most are restricted to a particular geographic region or habitat, such as soil, sediment, or water. To understand global trends, however, studies that analyze diversity across habitats and geographic regions are needed. Combining data from independent projects is oftentimes infeasible because community variation can be caused simply by differences in methodology. The Earth Microbiome Project (EMP) comprises 27,751 samples from 97 studies from a wide range of habitats and geographic regions that are processed in the exact same way. Although there are limitations to PCR-based sequencing surveys, this dataset is unique for its size in using standardize methods for all samples. Thus, this dataset provides an opportunity for a rigorous comparison of bacterial diversity across many parts of the globe.

3A recent overview from the EMP noted, as has previously been observed, that communities of free-living bacteria are more diverse than host-associated bacteria. For example, soil and sediment samples have higher alpha-diversity than animal gut or skin microbiomes. We expand on this initial alpha-diversity analysis by additionally evaluating beta- and gamma-diversity across spatial scales while mitigating for unevenly spaced samples (Fig 1). We also take the opportunity to use the large size of the EMP dataset to test whether different diversity metrics and taxa definitions influence our understanding of microbial diversity patterns.

4We specifically ask: which habitats support the highest levels of bacterial diversity? We consider three interrelated aspects of biodiversity: alpha-, beta-, and gamma-diversity. We measure alpha-diversity as the observed richness (number of taxa) or evenness (the relative abundances of those taxa) of an average sample within a habitat type. We quantify beta-diversity as the variability in community composition (the identity of taxa observed) among samples within a habitat. Finally, we calculate gamma-diversity as the total observed richness of all samples within in a habitat.

5We test several predictions about relative, not absolute, diversity patterns because, even in this large dataset, bacterial diversity remains undersampled. First, we predict that sediment and soil support the highest alpha-diversity within a single sample. These habitats are known to have relatively high bacterial diversity, although their relative rankings have not yet reached a consensus. Second, we expect that soil, sediment, inland water, and biofilm/mat habitats will exhibit high beta-diversity. These habitats are spatially separated with less dispersal or mixing than air or marine water. Finally, we predict that soils and sediments will exhibit high gamma-diversity as they are expected to have both high alpha- and beta-diversity.

6Within the soil habitat, we hypothesize that soils from biomes higher in plant diversity and productivity (e.g., forests and grasslands) support higher alpha-diversity than soils from biomes with low diversity and productivity (e.g., tundra and deserts). Of course, these biomes do not directly influence diversity, but they are defined based on abiotic factors, such as temperature or precipitation, that do influence diversity. Further, we expect that agricultural soils will exhibit lower beta-diversity than other biomes as common practices (pesticides, tilling, and fertilizer use) and the low diversity of crop plants influences community composition. We also compare the relationship between diversity and biomes to those between diversity, and pH or temperature to assess whether plant diversity or abiotic conditions more strongly influence bacterial diversity. We expect bacterial richness to peak at neutral pH and moderate temperatures, and abiotic factors to be a stronger influence on diversity than plant biomes. Overall, the aim of this study was to compare bacterial diversity trends across habitats. We show that the most diverse habitat depends on the type of diversity (alpha-, beta-, or gamma-diversity).

7Bacterial 16S rDNA (V4 region) sequence data and associated metadata (e.g., sample location, sample type, date of sampling) were downloaded from the Earth Microbiome Project (EMP) on September 1, 2016. Sample processing, sequencing, and core amplicon data analysis were performed by the Earth Microbiome Project (www.earthmicrobiome.org), and all amplicon sequence data and metadata have been made public through the data portal (qiita.microbio.me/emp). Data available from: https://doi.org/10.1038/nature24621. We used the EMP closed-reference (Greengenes 13.8) OTU dataset classified at 97% sequence similarity to reduce computational time (instead of the open-reference dataset). The dataset contains 27,751 samples, with a median depth of 54,091 sequences per sample. We excluded archaea from the analysis because, relative to bacteria, they make up a small portion of any given community (median = 0.018% of sequences).

8We used the EMP Ontogeny metadata to classify the habitat and, for soil, biome of each sample based on the EMP metadata (Fig 2). When the existing metadata were unclear, we used the latitude and longitude coordinates to assess the environmental context. Samples with insufficient data about their location or habitat were removed from the analysis. Host-associated samples were also removed. Further, we only retained samples that could be classified into one of the following habitats: soil, sediment, marine water, inland water (e.g., rivers and lakes), air, and biofilms/mats. These habitat types were chosen because they represent a wide range of environmental conditions and are well sampled within the EMP dataset. Within the soil habitat, we further classified samples into forest, hot desert, cold desert, grassland, shrubland, tundra, and agricultural soil. We also classified inland water and sediment samples as saline and non-saline. See supplemental materials for descriptions of sample locations (S1 Appendix). After removing samples with less than 15,000 sequences (rarefaction depth in this study), 11,680 free-living (non-host system) samples remained.

9To account for differences in sequencing depth, the samples were rarefied to 15,000 sequences with 1,000 resamplings in QIIME. This rarefaction depth provided a high sequence count per sample while minimizing sample loss to 4.74% of samples. All samples with less than 15,000 sequences were removed leaving 11,680 free-living (non-host system) samples. For each resampling, we calculated 24 alpha-diversity metrics on the rarified OTU table in QIIME (Fig 3A). These metrics characterized the community in five general ways: observed richness, estimated richness, evenness/dominance, phylogenetic diversity, and coverage of sampling. We used all 24 metrics throughout the alpha-diversity analysis to ensure that our final conclusions were not dependent on the type of metric. We calculated the median value of each metric across the 1,000 replicates.

10To minimize the effect of unevenly spaced samples, we averaged the alpha-diversity of the samples within a single geocluster. Many of the samples are highly clumped such that some geographic regions contribute unequally to the habitat’s diversity. Geoclusters (n = 172) were formed by clustering samples of the same habitat type located within 110 km of each other (distance of 1° latitude at the equator) using hclust and cutree from package ‘stats’ and rdist.earth from package ‘fields’ in R. While a smaller clustering distance would have yielded a higher geocluster sample size, this conservative distance allowed us to be more confident that our results reflected ecological processes, rather than sampling locations. We calculated the median of each diversity metric for the samples within each geocluster of the same habitat type. The averaged alpha-diversities were then cube root transformed to achieve normality and homoscedasticity. Finally, we tested for significant differences in alpha-diversity among habitats by performing a one-way ANOVA and Tukey’s HSD in R.

11We tested whether the alpha-diversity results depended on the diversity metric by running a correlation with every pairwise combination of diversity metrics using all 11,680 samples. Likewise, we tested whether our results depended on the resolution of OTU clustering by comparing the 97% similarity OTU table with the single-nucleotide resolution ‘sub-OTUs’ dataset, Deblur, produced by the EMP. We rarefied the Deblur dataset to 15,000 sequences per samples (1,000 times), calculated Exact Sequence Variance (ESV) richness per sample, and calculated the mean richness across the 1,000 replicates. We then ran a correlation between the OTU richness and ESV richness for every sample present in both datasets (n = 11,137).

12To further explore what factors might be driving alpha-diversity, we compared bacterial richness with pH and temperature at each sample site. We chose pH and temperature because these were the most widely included in the EMP dataset. For soil samples, however, the temperature data were often missing from the EMP dataset. Thus, for the analysis with just soil samples, we used temperature data from WorldClim. Metadata (pH and temperature) from the EMP were taken at the site at the time of sample collection. Data from WorldClim, a publicly-available dataset, included mean annual temperature averaged from 1970–2000 with spatial resolution of 10 minutes. Data is available from WorldClim Version2: http://doi.org/10.1002/joc.5086. We assigned external temperature data to soil samples using latitude and longitude with extract from package ‘raster’ in R. Temperature and pH were correlated with OTU richness using a second-degree polynomial in R. We tested whether temperature and pH differed among habitats and biomes using an ANOVA in R.

13To reduce computation time, we used a subset (150 rarefied tables) of the 1000 rarefied OTU tables generated during the alpha-diversity analysis to analyze beta-diversity. The OTU tables were first square root transformed to increase weight given to the rare taxa that make up the majority of microbial communities. For each of the 150 rarefied, square root transformed OTU tables, we calculated the median abundance for each taxon across all the samples within a single geocluster. We then calculated a Bray-Curtis dissimilarity matrix for each of the 150 OTU-by-geocluster tables in QIIME. Finally, we calculated the median of the 150 dissimilarity matrices to yield one median Bray-Curtis dissimilarity matrix.

14To visualize compositional differences among habitats, we used NMDS in PRIMER6. We also tested community composition differences among habitats using PERMANOVA in PERMANOVA+. Because we averaged OTU abundances for all samples of the same habitat type located within the same geographic area (geocluster), beta-diversity provides an approximation of the amount of community variation from location to location within one habitat (as opposed to variation from sample to sample within one location).

15To compare beta-diversity across habitats, we analyzed the variance within each habitat using the function PERMDISP in PERMANOVA+. To determine if unequal sampling among habitats biased these results, we re-calculated the Bray-Curtis values based on a selection of only 20 geoclusters for each habitat from a rarified OTU-by-geocluster table. We chose 20 geoclusters because that depth included five of the six habitats (excluding air) but avoided the biases expected with sample sizes less than ten. We repeated these subsamplings 100 times and tested for differences in beta-diversity among habitats using betadisper, the PERMDISP test implemented in the R package ‘vegan’. We compared the relative rankings of these rarefied beta-diversity results to the unrarefied results to determine if rarefaction changed the relationships of variance among habitat groups. Specifically, we considered the rarefied results to match the unrarefied results if 95–100 subsampled tests were significant and showed the same beta-diversity rankings (based on mean distance to centroid) as the unrarefied test.

16We tested whether the beta-diversity results depended on the diversity metric by running a correlation with every pairwise combination of nine diversity metrics. For each of the 150 OTU-by-geocluster tables, we calculated nine beta-diversity metrics using vegdist from package ‘vegan’ in R. We then took the mean matrix (of the 150 matrices) for each diversity metric. We performed a Spearman’s mantel test for every pairwise comparison of the nine averaged beta-diversity matrices using mantel from package ‘stats’. To compare Raup-Crick to the other beta-diversity metrics, we calculated the dissimilarity within, but not among, habitats. Raup-Crick is not an appropriate metric when communities do not share the same species pool. To calculate Raup-Crick, we took the mean of the matrices computed with raupcrick(…, chase = TRUE) and raupcrick(…, chase = FALSE) from package ‘vegan’ in R to follow the method recommended by Chase et al.. We calculated one Raup-Crick matrix for each habitat for each of the 150 OTU-by-geocluster tables, took the mean for each habitat across the 150 matrices, and then calculated the mean and SE distance from centroid in PRIMER6. We used the mean and SE to compare the trends in dissimilarity to those generated with the Bray-Curtis metric.

17To assess gamma-diversity by habitat, we plotted an OTU accumulation curve for each habitat with specaccum from package ‘vegan’ in R using the 150 OTU-by-geocluster tables. The OTU accumulation curve displays the numbers of geoclusters sampled on the x-axis and observed OTU richness on the y-axis. This plot allowed us to compare cumulative diversity levels across multiple samples distributed across the world. We examined whether habitats likely exhibit different gamma-diversity levels by calculating error bars equal to 1.96 times the standard deviation.

18Out of the six habitats compared, soils contained the highest observed richness (i.e., number of observed taxa rarefied at 15,000 sequences) for a single sample, with a median of 1,842 taxa (97% OTUs) per sample given this depth of sequencing (one-way ANOVA: F = 39.13, P < 0.001, r2 = 0.541; Fig 2A). Sediments were the second most diverse habitat with an average of 1,137 taxa. Marine water, air, inland water, and biofilms/mats had a significantly lower richness (averaging 571, 500, 478, and 342 taxa, respectively) than soils and sediments (P < 0.001; Fig 2A) but could not be distinguished from one another by richness.

19Because salinity influences bacterial community composition, we further tested whether taxon richness varied between non-saline and saline habitats. We found that salinity had no impact on alpha-diversity for sediments (one-way ANOVA: F = 0.433, P = 0.516; S1 Fig) or inland water (F = 0.093, P = 0.763; S1 Fig).

20Within the soil habitat, we further compared alpha-diversity among seven biomes (agricultural, grassland, shrubland, forest, hot desert, cold desert, and tundra soil). Within soil samples, richness differed significantly among biomes. Agricultural soils supported the highest richness in a sample, along with hot desert, grassland, and shrubland biomes. Forest soils were less diverse than agricultural soil, and tundra and cold deserts supported the lowest richness (one-way ANOVA: F = 42.62, P < 0.001, r2 = 0.642; Fig 2B). Notably, the cold desert biome was only represented by two geoclusters (averaging across 117 samples); thus, more data are needed to assess that particular biome’s diversity.

21The above results were robust to the alpha-diversity metric used. On a sample by sample basis, 24 alpha-diversity indices, including observed richness, were all correlated with each other (r2 = 0.09–1.00, P < 0.0001) with a mean r2 of 0.63 (Fig 3A). The metrics grouped into two main clusters. One cluster encompassed the richness/coverage metrics such as OTU richness, Faith’s Phylogenetic Diversity, and Chao1 (r2 = 0.86–1.00, mean r2 = 0.97). The other cluster included the evenness/dominance metrics such as Simpson’s and McIntosh dominance index (r2 = 0.31–1.00, mean r2 = 0.85). Further, each metric ranked the habitats from highest to lowest alpha-diversity in the same way, with the exception of air. Air communities were more even than other habitats, given their relative richness level (Fig 3B). Excluding air samples, OTU richness (ANCOVA: F = 54.68, P < 0.0001, r2 = 0.229), but not habitat (F = 2.14, P = 0.0775) was a predictor of evenness. When air samples were included, both OTU richness (F = 57.09, P < 0.0001, r2 = 0.173) and habitat (F = 7.0989, P < 0.0001, r2 = 0.107) were significant predictors of evenness.

22The alpha-diversity patterns were also robust to the taxonomic binning method (Fig 4). We compared the 97% OTU and Exact Sequence Variant (ESV) Deblur datasets provided by the EMP. Not only were OTU richness and ESV richness strongly correlated (r2 = 0. 933, P < 0.0001), but, on a sample by sample basis, they were also nearly identical (slope = 0.953, P < 0.0001). However, the relationship between OTU richness and ESV richness varied among habitats (ANOVA: F = 786.5, P < 0.0001, r2 = 0.331). In particular, non-saline sediments and inland water demonstrated a higher ESV:OTU richness ratio than other habitats (Fig 4).

23Taxon richness displayed a weak hump-shaped relationship with pH and a peak in diversity at a neutral pH (non-linear regression: P < 0.0001, r2 = 0.047; S2A Fig). In contrast, taxon richness only weakly correlated with temperature, and this relationship was driven by low-diversity biofilm/mat samples sampled from high temperatures (P < 0.0001, r2 = 0.036; S2C Fig). Overall, the bacterial alpha-diversity patterns across all habitats were not obviously related to pH or temperature. Most of the samples were, on average, at a neutral pH, with the soil samples more acidic and the biofilm/mat samples more basic (S2B Fig). Despite this, temperature (ANOVA: F = 272.3, P < 0.0001, r2 = 0.198) and pH (F = 245.2, P < 0.0001, r2 = 0.215) differed significantly among habitats (S2D Fig).

24Similar to the pattern observed across all habitats, richness within just the soil samples also peaked at a neutral pH (non-linear regression: P = 0.001, r2 = 0.121; S3A Fig). In contrast, richness in soils also peaked at a temperature around 10°C (P < 0.0001, r2 = 0.288; S3C Fig). Both temperature (ANOVA: F = 640.2, P < 0.0001, r2 = 0.153) and pH (F = 199.8, P < 0.0001, r2 = 0.594) differed among soil samples by biome (S3B and S3D Fig). Soils from both hot and cold deserts tended to be basic while agricultural fields, forests, and tundra were acidic. Tundra and cold deserts were the coldest biomes, and shrubland, agriculture, and hot deserts were among the hottest biomes.

25Sediment, biofilm/mat, and inland water habitats displayed the highest beta-diversity among geographic locations or geoclusters (not within a single sample), whereas soil, air, and marine water exhibited 17% lower beta-diversity (PERMDISP: F = 10.7, P = 0.001; Fig 2C). To test that these patterns were not influenced by unequal sampling (number of geoclusters) of the habitats, we subsampled the habitats (to 20 geoclusters per habitat) and retested the patterns. All 100 subsamplings produced the same beta-diversity rankings, and all models were significant, indicating that unequal sampling did not influence within-habitat beta-diversity. Within the soil habitat, beta-diversity did not differ by biome (P = 0.526; Fig 2D).

26Overall, bacterial community composition differed significantly by habitat (PERMANOVA: P = 0.001, Pseudo-F = 9.8601, r2 = 0.210; Fig 2E) and by biome for soils (P = 0.001, Pseudo-F = 4.221, r2 = 0.227; Fig 2F). Because salinity influenced the community composition for both sediments (P = 0.002, Pseudo-F = 2.0168, r2 = 0.075) and inland water (P = 0.02, Pseudo-F = 1.7702, r2 = 0.085), we tested whether salinity likewise influenced beta-diversity within these habitats. Beta-diversity did not differ between saline and non-saline samples within sediments (PERMDISP: P = 0.21, F = 2.4404) or inland water (P = 0.849, F = 0.27406; S4 Fig).

27The above results did not depend on the beta-diversity metric used. On a sample by sample basis, nine beta-diversity indices, including Bray-Curtis, were correlated with each other (r = 0.435–1.00, P = 0.001, mean r = 0.88; S5A Fig). Raup-Crick has been suggested as a more appropriate metric when comparing groups with different alpha-diversity levels. Because Raup-Crick assumes that all communities are part of the same regional species pool, we calculated the mean and standard error within each habitat (excluding between habitat comparisons) and compared the trend to that generated by the Bray-Curtis metric. Both Bray-Curtis and Raup-Crick metrics showed the same trend of beta-diversity among habitats (S5B Fig).

28Considering the accumulation of taxon richness across geoclusters, sediments exhibited the highest gamma-diversity of any habitat, followed by soils and inland water (Fig 5). The sediment rarefaction curve showed little sign of flattening out, indicating that most taxa are yet to be sampled. In contrast, the soil curve noticeably leveled off, even at a similar level of sampling. The gamma-diversity of marine water, biofilms/mats, and air were not statistically distinguishable from one another but, as a group, exhibited lower gamma-diversity than inland water, soils, and sediments.

29Here, we tested which habitat contains the most bacterial taxa within a single sample (alpha-diversity), which exhibits the most variation among samples (beta-diversity), and which contains the most taxa across all samples (gamma-diversity). We show that a single sample of soil on average contained higher bacterial alpha-diversity than any other habitat, including sediment (Fig 2A). However, sediment had higher gamma-diversity, with much of its diversity yet to be sampled (Fig 5). Within soils, we found that agricultural soils had among the highest richness and exhibited just as much compositional variation (beta-diversity) as other biomes.

30Although both sediments and soil were previously known to be highly diverse microbial habitats, previous studies demonstrated conflicting results about their relative ranking. Using 11,680 samples and minimizing geographic biases, this analysis suggests that soil contains higher alpha-diversity than sediments. In contrast, marine water, inland water, air, and biofilms/mats contain the lowest alpha-diversity (Fig 2A).

31Within-habitat heterogeneity is a known driver of plant and animal diversity and has been correlated with microbial communities as well. Here, we show that the alpha-diversity patterns are consistent with the idea that habitat heterogeneity may drive bacterial diversity at a single sample. The highly mixed water and air environments harbor lower diversity, consistent with previous smaller-scale studies. While both sediments and soil are not as well mixed, sediments contain higher water content than soils. Water content increases connectivity and thus reduces environmental heterogeneity and promotes dispersal, both of which can result in lower diversity. At the same time, biofilms and mats, despite being spatially structured, also displayed low alpha-diversity. However, the biofilm/mat samples from the EMP dataset encompassed samples with the highest pH and temperature (S2 Fig). We therefore speculate that these abiotic extremes contribute to low alpha-diversity. In fact, diversity in many habitats is lowest at extreme temperatures. Yet ultimately, little is known about the environmental conditions, and their heterogeneity, at the spatial scale that matters for microorganisms. To test the importance of within-sample heterogeneity on microbial diversity directly, finer-scale data are needed.

32Our analysis is the first to quantify bacterial beta-diversity among habitats across many parts of the globe. While soils contained the highest alpha-diversity within a single sample, sediments displayed higher beta-diversity among geoclustered samples within a habitat (Fig 2C). Sediment beta-diversity was also similar to that of inland water and biofilms/mats. Additional environmental data associated with the individual samples would be needed to distinguish whether these beta-diversity patterns might be driven by dispersal limitation or spatial variation in environmental conditions.

33Given their high alpha- and beta-diversity, it is not surprising that sediments are also estimated to contain the highest gamma-diversity (Fig 5). While extracellular DNA (eDNA) may be particularly prevalent in ocean sediments, evidence thus far suggests that eDNA has minimal effect on sediment diversity estimates from sequencing surveys. The taxa accumulation curves also suggest that, while we may have observed most bacterial taxa in soil (at least from highly sampled continents), there is much more diversity to discover in sediments. Similarly, other than soil, the accumulation curves suggest that air, water habitats, and biofilms/mats remain undersampled as well. Of course, the number and localities of samples available will influence the diversity estimates. Thus, a limitation to these conclusions is that the samples are highly concentrated in North America and Europe (Fig 1), and continued sampling is needed to test the robustness of these diversity patterns.

34Because plant diversity and productivity are shown to impact microbial communities, we further characterized alpha- and beta-diversity trends among biomes from which the soil samples were collected. Agricultural soils contained among the highest alpha-diversity, as previously noted in smaller scale studies. Indeed, some agricultural practices, such as application of manure, are known to increase bacterial diversity. Even more notable, however, is that agricultural soils encompassed similar levels of beta-diversity to those of other biomes. While some agricultural practices have been shown to homogenize communities within a single field, not all practices have a homogenizing effect. Further, the diversity of agricultural practices around the world seems to select for as much variation in bacterial composition (beta-diversity) as different types forests or deserts.

35Contrary to our hypothesis, biomes with higher plant diversity or productivity, such as forest or shrubland soils, were no more diverse within a sample than other biomes (Fig 2B). These results support previous findings that grasslands contain more bacterial diversity than forests and, overall, plant and soil diversity are uncoupled. We therefore propose that abiotic factors may be more important for soil bacterial alpha-diversity than plant biomes. Biomes differed significantly in pH and temperature, and soil alpha-diversity was strongly correlated with both factors (S3 Fig). These results are consistent with previous studies that find bacterial richness in soils peaks at a neutral pH and at mid-temperatures (around 10°C). Of course, other unmeasured environmental factors and/or ecological interactions are likely influencing soil diversity.

36Finally, these diversity patterns appear to be robust to two key methodological issues. First, diversity trends did not depend on the particular alpha- or beta-diversity metrics used (Figs 3A and S5A). Air, as the only exception, was more even than expected, given its richness (Fig 3B). We speculate that the movement of air contributes to its evenness as air likely picks up a sampling of bacteria from many different habitats. Second, the results were robust to the degree of clustering of the amplicon sequences (Fig 4). Both the 97% OTU and the ESV datasets yielded the same alpha-diversity trends, as previously noted in a smaller scale study. Most of those outliers in this analysis came from non-saline sediment or inland water samples, which had higher ESV richness than 97% OTU richness. While this pattern could suggest higher finer-scale diversity within these habitats, we caution that these samples originated from only four geoclusters. Overall, while ESVs can be useful for resolving finer diversity among specific taxonomic groups, broad-scale alpha-diversity patterns do not seem to be altered by these particular operational definitions.

37With the largest dataset created with consistent methodology and a geographically widespread sampling effort, we show that soils support the highest diversity within a single sample (alpha-diversity) and that sediments are more variable in composition among locations (beta-diversity) and likely support the most bacterial taxa at a larger spatial scale (gamma-diversity). Within soils, we find biome type impacts soil alpha-diversity but not beta-diversity. Many of these results appear consistent with the idea that spatial heterogeneity and dispersal limitation promote bacterial diversity. These baseline patterns set the stage for new research on the mechanisms driving the generation and maintenance of bacterial diversity.