Global Microbiome Diversity Scaling in Hot Springs With DAR (Diversity-Area Relationship) Profiles

The spatial distribution of biodiversity (i.e., the biogeography) of the hot-spring microbiome is critical for understanding the microbial ecosystems in hot springs. We investigated the microbiome diversity scaling (changes) over space by analyzing the diversity-area relationship (DAR), which is an extension to classic SAR (species-area relationship) law in biogeography. We built DAR models for archaea and bacteria with 16S-rRNA sequencing datasets from 165 hot springs globally. From the DAR models, we sketch out the biogeographic maps of hot-spring microbiomes by constructing: (i) DAR profile—measuring the archaea or bacteria diversity scaling over space (areas); (ii) PDO (pair-wise diversity overlap or similarity) profile—estimating the PDO between two hot springs; (iii) MAD (maximal accrual diversity) profile—predicting the global MAD; (iv) LRD/LGD (ratio of local diversity to regional or global diversity) profile. We further investigated the differences between archaea and bacteria in their biogeographic maps. For example, the comparison of DAR-profile maps revealed that the archaea diversity is more heterogeneous (i.e., more diverse) or scaling faster than the bacterial diversity does in terms of species numbers (species richness), but is less heterogeneous (i.e., less diverse) or scaling slower than bacteria when the diversity (Hill numbers) were weighted in favor of more abundant dominant species. When the diversity is weighted equally in terms of species abundances, archaea, and bacteria are equally heterogeneous over space or scaling at the same rate. Finally, unified DAR models (maps) were built with the combined datasets of archaea and bacteria.

Introduction

Hot springs are one of the extreme environments on the earth planet. Hot spring microbiomes play a critical role in shaping the geothermal ecosystems. The structures and functions of microbial communities inhabiting hot springs have their somewhat unique characteristics compared with non-geothermal environment microbiomes (Inskeep et al., 2013a,b; Song et al., 2013; Masaki et al., 2016; Poddar and Das, 2018), such as soil microbiome (Fierer et al., 2012; Hartmann et al., 2014), marine microbiome (Gajigan et al., 2018), and human microbiome (Huttenhower et al., 2012). Hot springs are often abundant in thermophilic, hyperthermophilic, and thermoresistant bacterial and archaeal taxa (Urbieta et al., 2014a,b). For example, the hot springs located in the Chilas and Hunza areas of Pakistan with 90–95°C temperature harbor abundant phylum Thermotogae (Amin et al., 2017). Many reports have suggested that hot spring microbial communities are extremely heterogeneous and are often dominated by the thermophilic bacterium (e.g., Cole et al., 2013; Sharp et al., 2014; Masaki et al., 2016; Poddar and Das, 2018). Even though some recent studies showed that microbial eukaryotes, especially microbes from the phyla of Ascomycota and Basidiomycota, may also be important components of hot springs microbiomes (de Oliveira et al., 2015; Salano et al., 2017; Liu et al., 2018; Oliverio et al., 2018), the microbial communities routinely consist of Bacteria and Archaea (Meyer-Dombard and Amend, 2014; Hedlund et al., 2015; Merkel et al., 2017).

The stability of hot spring environments is routinely determined by the steady state of their microbial diversity in a specific environment, where water temperature, pH, and chemical composition are often the most important factors to influence the diversity (Mohanrao et al., 2016; Amin et al., 2017; Chan et al., 2017; Ghilamicael et al., 2017; Poddar and Das, 2018). In general, there is an inversely proportional relationship between microbial diversity and temperature of hot spring (Cole et al., 2013; De León et al., 2013; Amin et al., 2017; Chan et al., 2017). In addition to temperature, water pH is another primary environmental factor directly influencing microbial diversity in hot springs (Inskeep et al., 2013a,b; Xie et al., 2015). The water pH is determined by the chemical composition in hot springs. The role of chemical composition in designing the structure of geothermal microbial communities should not be underestimated, which sometimes play together with the pH (Jiang et al., 2015; Geesey et al., 2016). In spite of the extensive studies on the microbial diversities in hot springs, which suggest that the biodiversity of microbial communities vary with physicochemical conditions and biogeographical location of inhabiting hot springs, to the best of our knowledge, the hot-spring microbiome diversity scaling on regional or global scales from a biogeography perspective has not been investigated yet.

To investigate the microbiome diversity scaling on regional/global scale, one of the most powerful theoretical tools is the classic species-area relationship (SAR) power law, which achieved a rare “law” status in ecology and biogeography. The SAR is often described with a power function S = cA^z, where S is the number of species accumulated from a region of size A, and z is termed species (number) scaling parameter. The study of SAR can be traced back to the nineteenth century (Watson, 1835; Arrhenius, 1921; Preston, 1960, 1962) and the relationship was said to inspire MacArthur and Wilson (1967) to establish their island biogeography theory, which helped to shift the focus of ecology from population level to community level.

Theoretically, a series of extensions of the classic SAR to general diversity-area relationship (DAR) were introduced by the author's group recently (Ma, 2018a,b). The extensions were justified to remedy a limitation with classic SAR, where the biodiversity is measured with the number of species or the so-termed species richness. Species richness can be a rather meaningful measure for biodiversity in the case of large plants and animals, but in many other cases (especially for microbes), it is a poor measure of biodiversity because it ignores the differences in species abundances. For example, 1,000 of panda and one billion of panda will weigh in the same with species richness; but if the latter number were the case, panda would not have been on the endangered species list. The DAR extension was facilitated by the adoption of the Hill numbers (Hill, 1973; Chao et al., 2012, 2014) as general diversity measures, which weight diversity differently depending on the so-termed diversity order. In terms of the Hill numbers, biodiversity can be measured by the so-termed diversity profile (Chao et al., 2012, 2014), which calculates a series of Hill numbers, weighted differently by the species abundance distribution (SAD). Therefore, Hill numbers are now well-recognized as the most appropriate measures for alpha-diversity and its multiplicative partition of beta-diversity is also considered with advantages over other beta-diversity measures. With the new DAR approach, four sets of new tools: DAR profile, PDO (pair-wise diversity overlap) profile, and MAD (maximal accrual diversity) profile, and LRR/LGR [local to regional (global) diversity ratio], can be established with the parameters from DAR modeling. These profiles, together with DAR models offer powerful tools, not only for quantifying the regional/global scaling of biodiversity, but also for sketching out the biogeography maps of biodiversity distribution (Ma, 2018a,b). In the present study, we apply the DAR approach to analyzing the global biodiversity scaling of hot-spring microbiome by reanalyzing the 16S-rRNA marker gene abundance datasets of 165 hot springs on a global scale, previously collected and published by Sharp et al. (2014). We further sketch out and compare the biogeography “maps” of archaea and bacteria, and highlight their differences in biodiversity distribution.

The four profiles, DAR (diversity-area relationship), PDO (pair-wise diversity overlap), MAD (maximal accrual diversity) profile, LRD (local to regional diversity ratio), we build for archaea, bacteria, and their combined assemblages offer tools to sketch out the biogeography maps with different themes. While, the map theme profiled by DAR is the diversity scaling (difference or heterogeneity) over space, the theme profiled by PDO is the similarity of diversity over space. While the map theme profiled by MAD is the theoretically maximal accrual diversity (essentially the maximal gamma diversity), the theme profiled by LRD/LGD is then the local vs. regional diversity comparison, which answers a simple question: how much, on average, a local sample can represent the regional or global diversity.

Materials and Methods

The Hot Spring Microbiome Datasets

The datasets of 165 hot-spring microbiome were originally collected and reported by Sharp et al. (2014). Their 16S-rRNA OTU (operational taxonomic unit) tables were generated from 165 microbiome samples taken from sediment, soil, and mat in Western Canada and Taupo Volcanic Zone, New Zealand (Sharp et al., 2014). A total of 1,162,553 high quality sequences were obtained from the 165 samples with 634–15,283 sequences per sample. There were 61,910 OTUs, including 7,964 archaea and 53,946 bacteria, when those sequences were clustered at the 97% identity threshold. Further information on the datasets is referred to Sharp et al. (2014).

Definitions and Computational Procedures

Three steps are involved in building DAR models with microbiome datasets (see Figure 1): (i) bioinformatics analysis of 16S-rRNA data to get OTU tables (e.g., Schloss et al., 2009; Caporaso et al., 2010; Sinha et al., 2015). The microbiome quality control project: baseline study design and future directions. Genome Biology. Vol. 16: 276, https://doi.org/10.1186/s13059-015-0841-8; (ii) computing species or OTU diversities measured with the Hill numbers (Chao et al., 2012, 2014; Ma, 2017); (iii) constructing the DAR models (Ma, 2018a,b).

FIGURE 1

Figure 1. A diagram showing the three steps involved in DAR (diversity-area relationship) analysis: major software packages (MOTHUR, QIIME, QIIME2), data format (OTU table), definitions (Hill numbers), models (PL and PLEC), and concepts (profiles) were illustrated.

Diversity Measured in Hill Numbers

The Hill numbers (Hill, 1973; Jost, 2007; Chao et al., 2012) are considered as the most appropriate measure for alpha diversity,

$\begin{array}{ll}{}^{q}D={\left({\displaystyle \sum _{i=1}^S}{p}_i^q\right)}^{1/(1-q)}& (1)\end{array}$

where S is the number of species, p_i is the relative abundance of species i, q is the order number of diversity. The Hill number is not defined when q = 1, but its limit as q approaches to 1 exists in the following form:

$\begin{array}{ll}{}^{1}D={{\displaystyle \underset{q\to 1}{lim}}}^{q}D=exp\left(-{\displaystyle \sum _{i=1}^s}{p}_{i}log(p_1)\right)& (2)\end{array}$

The parameter q determines the sensitivity of the Hill number to the relative frequencies of species abundances. If q = 0, the species abundances do not weigh at all and ⁰D = S, which is simply the species richness. When q = 1, ¹D equal the exponential of Shannon entropy, and is interpreted as the number of typical or common species in the community because ¹D is weighted proportionally by species abundances. When q = 2, ²D equal the reciprocal of Simpson index, i.e., the number of dominant or very abundant species in the community (Chao et al., 2012) because ²D is weighted in favor of more abundant species. The general interpretation of ^qD (diversity of order q) is that the community has a diversity of order q, which is equivalent to the diversity of a community with ^qD = x equally abundant species. The so-termed diversity profile refers to the Hill numbers at different diversity order q (Jost, 2007; Chao et al., 2012, 2014).

The DAR Models and DAR Profile

Ma (2018a) extended SAR (species area relationship) to general DAR (diversity area relationship), in which diversity is measured with Hill numbers. The first DAR model, which borrowed the same power law (PL) function from the classic SAR, is:

$\begin{array}{ll}{}^{q}D=c{A}^z& (3)\end{array}$

where ^qD is diversity measured in the q-th order Hill numbers, A is area, and c and z are parameters.

A second DAR model is the power law with exponential cutoff (PLEC) model, which was originally introduced to SAR modeling by Plotkin et al. (2000), Ulrich and Buszko (2003), and Tjørve (2009) is:

$\begin{array}{ll}{}^{q}D=c{A}^{z}exp(dA),& (4)\end{array}$

where d is a third parameter and is usually negative in the DAR models, and exp(dA) is the exponential decay term, which eventually overwhelms the power law behavior at very large value of A.

The following log-linear transformed equations can be used to estimate the parameters of the DAR models:

$\begin{array}{ll}\text{ln}(D)=\text{ln}(c)+z\text{ ln}(A)& (5)\end{array}$

$\begin{array}{ll}\text{ln}(D)=\text{ln}(c)+z\text{ ln}(A)+dA& (6)\end{array}$

The linear correlation coefficient (R) and p-value are used to judge the goodness of the model fitting.

Ma (2018a) defined the relationship between DAR-PL (power law) model parameter (z) and diversity order (q), or z-q trend, as the DAR profile, which comprehensively describes the change of diversity scaling parameter (z) with the diversity order (q).

Predicting MAD (Maximal Accrual Diversity) With PLEC-DAR Models

Ma (2018a) derived the maximal accrual diversity (MAD) in a cohort or population based on PLEC model [Equations (4) and (6)] as follows:

$\begin{array}{ll}Max{(}^{q}D){=}^q{D}_{max}=c{(-\frac{z}{d})}^z\text{exp}(-z)=c{A}_{max}^z\text{exp}(-z)& (7)\end{array}$

where A_max is the number of areas accrued to reach the maximum and is equal to:

$\begin{array}{ll}{A}_{max}=-z/d& (8)\end{array}$

and all parameters are the same as in Equations (4) and (6).

Similar to the previous definition for DAR profile (z-q pattern), Ma (2018a) defined the MAD profile (D_max-q pattern) as a series of D_max values corresponding to different diversity order (q).

Pair-Wise Diversity Overlap (PDO) Profile

The pair-wise diversity overlap (g) of two bordering areas of the same size (i.e., the proportion of the new diversity in the second area) is:

$\begin{array}{ll}g=2-2^z& (9)\end{array}$

where z is the scaling parameter of DAR-PL model [Equations (3) and (5)]. If z = 1, then g = 0 and there is no overlap (similarity); and if z = 0, then g = 1 and there is a total overlap. In reality, g should be between 0 and 1. Since g is between 0 and 1, one may even use percentage notation to measure PDO.

Similar to previous definitions for DAR profile (z-q pattern) and MAD profile (D_max-q pattern), Ma (2018a) defined the PDO profile (g-q pattern) as a series of PDO values corresponding to different diversity order (q).

The Ratio of Local Diversity to Regional (or Global) Accrual Diversity

Ma and Li (2019) defined the LRD (or LGD) as the ratio of local diversity of an averaged area to the regional diversity accrued or the global MAD (maximal accrual diversity). The dividend (local diversity) is ideally estimated with the parameter c of the DAR-PL model, but can be approximated with the parameter c of the DAR-PLEC model. The divisor can be either regional accrual diversity (which can be estimated with PLEC model directly) or global maximal accrual diversity (which is simply the MAD or D_max). Hence, in general, two similar metrics can be defined, depending on the regional or global scale is adopted: one is the ratio of local diversity to regional diversity (LRD), and another is the ratio of local diversity to global MAD (LGD). The LRD (LGD) can be computed with the following formulae, respectively:

$\begin{array}{ll}LRD=c/D& (\mathrm{10a})\end{array}$

$\begin{array}{ll}LGD=c/D_{max}& (\mathrm{10b})\end{array}$

where D can D_max be computed with the PLEC model directly (Equations 4 and 7, respectively), c can be estimated or approximately with the PL or PLEC model. The LRD (LGD) at different diversity orders (q = 0, 1, 2, …) were defined as LRD (LGD) profile, or local to regional (local to global) diversity scaling profile (Ma and Li, 2019). It is essentially the ratio of alpha to gamma diversity.

Re-sampling Procedure to Enhance the Robustness of DAR Modeling

The accumulation order of areas in DAR modeling may influence the estimation of parameter c in fitting PL/PLEC models (Equations 3–6). When there is not a natural spatial sequence (or arrangement) among the communities sampled, or the arrangement information is not available, arbitrarily choosing an accumulation order (arrangement) can be problematic. To avoid the potential bias from an arbitrary order of the hot spring microbiome samples, we totally permutated the orders of all the community samples under investigation, and then randomly choose 100 orders of the communities generated from the permutation operation. In other words, rather than taking a single arbitrary order for accruing community samples in one-time fitting to the DAR model, we iteratively perform the DAR model-fitting 100 times with the 100 randomly chosen orders. Finally, the averages of the model parameters from the 100 times of DAR fittings are adopted as the model parameters of the DAR. In the case of this study, we do not have detailed information on the geographic locations of the hot-spring microbiome samples, the re-sampling scheme is adopted to remedy the deficiency.

Results and Discussion

DAR Analysis of the Archaea

Table 1 displays the parameters from fitting the DAR (diversity-area relationship) PL (power law) and PLEC (power law with exponential cutoff) models to the datasets of archaea in hot springs. The p-values in Table 1 show statistically significant fitting (in all models p < 0.001) of both PL and PLEC to the datasets. The PLEC model has an advantage of being able to estimate the MAD (maximal accrual diversity) or D_max, which essentially measures the accrued diversity in a population or cohort, with the so-termed MAD-profile, as explained previously. The other two parameters from PL model, scaling parameter (z) and pair-wise diversity overlap (PDO) parameter (g) define the DAR-profile and PDO profile, respectively. From Table 1, we summarized the following findings:

   (i) The archaea diversity scaling (changes across space) fitted to the PL model successfully, with p < 0.001, and succeeded in all 100 times of re-sampling. It also fitted to PLEC model successfully p < 0.001, but with slightly less robustness given that some failures out of the 100 times of re-sampling occurred, especially at the high diversity orders. This does not imply that PLEC model is not applicable to the archaea DAR at all; it only suggests that the order of accruing areas (hot spring samples) influences the fitting of the model.

 (ii) The DAR profile: the parameter z-q series from the PL model is: z-q = {0.867[q = 0], 0.575[q = 1], 0.512[q = 2], 0.492[q = 3]}, a monotonically decreasing trend with the increase of diversity order (q). This indicates that at lower diversity orders, the spatial heterogeneity is larger than that at higher orders. The highest heterogeneity is at the species richness level (q = 0). Since at higher diversity order, the dominant species were given more weights in computing the diversity (Hill numbers). This suggests that the hot spring microbiomes are more homogenous (less heterogeneous or less diverse) in terms of their dominant species. Since at lower diversity order, the rare species were given more weights in computing the diversity (Hill numbers). This suggests that the hot spring microbiomes are more heterogeneous (more diverse) in terms of rare species.

(iii) The PDO profile: the parameter g-q series from the PL model is: g-q = {0.170[q = 0], 0.505[q = 1], 0.568[q = 2], 0.588[q = 3]}, a monotonically increasing trend with the increase of diversity order (q). This indicates that at lower diversity orders, the pair-wise diversity overlap (similarity) is smaller than that at higher orders. The lowest diversity overlap (similarity) is at the species richness level (q = 0), and the highest similarity occurred at the diversity order q = 3. This finding further supports the finding revealed by the DAR-profile in the previous item (ii).

 (iv) The MAD profile: the parameter D_max-q series estimated from the PLEC model is: D_max-q = {8,397[q = 0], 248.3[q = 1], 68.8[q = 2], 45.0[q = 3]}, shows a monotonically decreasing trend with the increase of diversity order (q). This, of course, is determined by the nature of the diversity (Hill numbers). The D_max(q = 0) = 8,397 is simply the maximal accrual of species richness since the Hill numbers at q = 0 is the species richness. For example, at diversity order q = 1, the Hill numbers (our measure for the diversity) are equivalent to (function of) Shannon diversity index. This suggests that the maximal accrual diversity in terms of the Shannon index is 248.3. Similar interpretations can be made for q = 2, and 3. The maximal accrual species richness of 8,397 means that, it needs to accrue A_max = 262 hot spring sites (areas) globally to reach this theoretical asymptote of species richness in the hot spring microbiome.

  (v) There appears a trend of decreasing correlation coefficient (R) with increasing diversity order (q). This should be expected because with increasing q, the complexity associated with non-linearity in higher order entropy (i.e., Hill numbers) is increased. Consequently, goodness-of-fitting to the linear models (Equations 5, 6) is likely to decline.

TABLE 1

Table 1. The parameters of DAR (diversity-area relationship) for the Archea in the hot springs, computed with 100 times of re-sampling from the totally permutated hot springs.

DAR Analysis of the Bacteria

We did the same DAR analysis with bacteria dataset, and the results are exhibited in Table 2. We further performed Wilcox non-parametric significance test of the differences between Archaea and Bacteria in their DAR parameters, and it turned out that (i) regarding the PL-model, archaea-DAR and bacteria-DAR have significantly different DAR parameter values except for the diversity order q = 1. (ii) Regarding the PLEC model, archaea-DAR, and bacteria-DAR have significantly different DAR parameter values at the higher diversity orders (q = 2, 3), but no significant differences occurred at the lower diversity orders (q = 0, 1). These test results justify our attempt to separately build DAR models for archaea and bacteria. Since the format of Tables 1, 2 are exactly the same, our explanations for the bacteria-DAR is presented relatively brief intentionally. From Table 2, we summarized the following findings:

   (i) The bacteria diversity scaling (changes across space) fitted to both the PL and PLEC models successfully, with all p < 0.001. The goodness-of-fitting is equally well with that for the archaea DAR models in the previous sub-section.

 (ii) The DAR profile: the parameter z-q series from the PL model is: z-q = {0.830[q = 0], 0.616[q = 1], 0.575[q = 2], 0.544[q = 3]}, a monotonically decreasing trend with the increase of diversity order (q). This indicates that at lower diversity orders, the spatial heterogeneity is larger than that at higher orders. The highest heterogeneity is at the species richness level (q = 0). This pattern per se is the same as that of archaea-DAR.

If we further compare both the DAR profiles (see statistical tests in Table 3), we found that archaea has a larger diversity scaling parameter (z-values) at diversity order q = 0 (i.e., species richness, but smaller scaling parameter (z-values) at diversity order q = 2 or 3. At diversity order q = 1, which is equivalent to the diversity measured with Shannon entropy and weighs all species in proportion with their relative abundance levels, archaea and bacteria showed no significant difference in their scaling parameter (z-values). These findings indicate that archaea is more heterogeneous or scaling faster than bacteria does in terms of species numbers (species richness), but is less heterogeneous or scaling slower than bacteria when the diversity (Hill numbers) were weighted computationally in favor of more abundant dominant species. When the diversity (Hill numbers) is weighted equally in terms of species abundances, archaea, and bacteria are equally heterogeneous over space or scaling at the same rate.

The above finding also highlighted the necessity of using Hill numbers as general diversity measures (the diversity profile) over using a single ad-hoc diversity measure such as Shannon entropy or Simpson's index, because the latter may lead to inconsistent results or loss of information. This also shows the necessity of using DAR profile, a series of scaling parameter values (z) across different diversity orders (q), rather than using a single scaling parameter as in classic SAR (species-area relationship) analysis.

(iii) The PDO profile: the parameter g-q series from the PL model is: g-q = {0.217[q = 0], 0.456[q = 1], 0.502[q = 2], 0.535[q = 3]}, a monotonically increasing trend with the increase of diversity order (q). This pattern is the same as that of the archaea DAR PDO profile.

Similar to the previous DAR-profile comparison between archaea and bacteria, the archaea has a smaller PDO overlap (similarity) than the bacteria has at species richness level (q = 0), but has a larger PDO overlap (similarity) at the higher diversity order q = 2 or 3. At the diversity order q = 1, archaea and bacteria have the same level of diversity overlap (similarity) across space. The interpretation for this finding is exactly the same as that for the DAR-profile above.

 (iv) The MAD profile: the parameter D_max-q series estimated from the PLEC model is: D_max-q = {55,489[q = 0], 2,831.1[q = 1], 427.7[q = 2], 207.9 [q = 3]}, showing a monotonically decreasing with the increase of diversity order (q). This, of course, is determined by the nature of the diversity (Hill numbers). The D_max(q = 0) = 55,489 is simply the maximal accrual of species richness since the Hill numbers at q = 0 is the species richness. Similar interpretations can be made for q = 1, 2, and 3. The maximal accrual species richness of 55,489 means that, it needs to accrue A_max = 256 hot spring sites (areas) globally to reach this theoretical asymptote of species richness in the hot spring microbiome.

  (v) Similar to the previous DAR models for archaea, there appears a trend of decreasing correlation coefficient (R) with increasing diversity order (q). This should be expected as explained previously.

TABLE 2

Table 2. The parameters of DAR (diversity-area relationship) for the Bacteria in the hot springs, computed with 100 times of re-sampling from the totally permutated hot springs.

TABLE 3

Table 3. The p-values of Wilcox non-parametric significance test between the differences between Archaea and Bacteria in their DAR parameters.

While the above pattern for bacteria MAD-profile is similar to the pattern for archaea MAD profile, the vis-a-vis comparison of both the MAD profiles is a different story. Obviously, the values of bacteria D_max are far larger than the values of archaea D_max. Indeed, the difference is consistent with biological (ecological) reality that there are far more bacteria species than archaea species in hot springs. Unfortunately, unlike the cases of DAR and PDO profiles, we cannot perform the permutation (randomization) test for MAD-profile (D_max). This is because the D_max was computed based on the average parameter values form 100 times of re-sampling. We believe biological (ecological) observations justify our claims that MAD-profiles are also different between the archaea and bacteria.

DAR Analysis With the Combined Datasets of Archaea and Bacteria

Since there are significant differences between the archaea and bacteria in their DAR parameters, ideally, independent DAR models should be built for each kingdom. However, there is no doubt that they cohabitate (coexist) in the hot spring environment. Therefore, building unified DAR models (Table S7) for the combined archaea and bacteria is justified. Table S7 shows that the DAR models fitted to the combined datasets of archaea and bacteria equally well with those for the archaea or bacteria, independently. For practical purpose such as conservation planning, the unified models (Table S7) are obviously more convenient, but for theoretical (mechanistic) inquiries, the separately built DAR models previously (Tables 1, 2) should be more appropriate. Since the pattern of the unified DAR models are similar to the separately built ones, except some nuances, which make little differences for practical applications. As to the theoretical implications of those nuances, we recommend the use of those separately built DAR models directly. Therefore, we do not further compare the subtle differences between the unified and separate DAR models here.

The Ratio of Local Diversity to Regional (or Global) Accrual Diversity

The LRD (or LGD) is defined as the ratio of the local diversity of an averaged area to the regional diversity [or the global maximal accrual diversity (MAD)]. The dividend (local diversity) can be estimated with parameter c of the DAR-PL model, and the divisor can be either regional accrual diversity (which can be estimated with PLEC model directly) or global maximal accrual diversity (which is simply the MAD or D_max). Note that we defined LRD/LDG at different diversity orders (q = 0, 1, 2, 3, …) as LRD/LGD profile.

Here we only computed LGD, the global version of the ratio. Table 4 listed the LGD for the archaea, bacteria and their combined microbiome, at each diversity order (q). For example, at species richness level (q = 0), the LGD is between 1.24 and 1.57%, which suggests that, on average, a single (local) hot spring only hosts between 1.24 and 1.57% of the global scale diversity. At high diversity orders, the ratios increased (up to 9% approximately). Another interesting observation is that at the species richness (q = 0), the LGD for archaea is lower than that for bacteria. However, at higher diversity orders (q = 1, 2, 3), the trend is reversed.

TABLE 4

Table 4. The LGD (the ratio of local diversity to global maximal accrual diversity) profile for the archaea, bacteria, and combined communities in the hot springs.

Discussion

With the gold rush of microbial community ecology, thanks to the revolutionary metagenomic sequencing technology, the classic SAR has been called for new missions. Green et al. (2004) and Horner-Devin et al. (2004) published, in the same issue of the journal Nature, the first two studies on the SAR of microbes. The following year, two other important studies by Bell et al. (2005) and Smith et al. (2005) were published in two other leading journals, Science and PNAS, respectively. The SAR power law exponent (b) values from those studies were 0.074 (fungi), 0.019–0.040 (bacteria in marsh sediment), 0.26 (bacteria in tree holes), and 0.134 (phytoplankton). According to Green and Bohannan (2006) review, the reported SAR exponents in microbes were in the range between 0.019 and 0.470, but most values were below 0.2 (8 out of 11 studies). A major limitation of these pioneering studies on the testing of SAR with microbes is then low throughput of DNA sequencing technology in detecting bacteria, and consequently the diversity and SAR exponent were significantly underestimated. Even with the technology limitation, the reported exponent values have already indicated the applicability of SAR in microbes, and recent studies further confirmed the validity of microbial SAR (e.g., Noguez et al., 2005; Peay et al., 2007; Bell, 2010; Barreto et al., 2014; Pop Ristova et al., 2014; Ruff et al., 2015; Terrat et al., 2015; Várbíró et al., 2017). For example, nearly a decade after Green and Bohannan (2006) review, the range of exponent (z) of microbial SAR is nearly unchanged and most studies have still been limited to bacteria and archaea (Barreto et al., 2014).

While the classic SAR has been well-recognized as one of the most significant laws in ecology and biogeography, it is not without limitations. The recent extension from SAR to DAR by Ma (2018a,b) generalized the scaling law of biodiversity from species richness (the number of species) to general diversity measures (the Hill numbers). Furthermore, DAR profile, PDO profile, MAD profile, LRD/LGD profile based on DAR models can offer useful novel tools to sketch out the biogeography maps, which comprehensively characterize the biodiversity scaling over space and time.

Despite the large number studies of microbial SAR in various environments, to the best of our knowledge, the SAR of the hot-spring microbiome has not been reported in existing literature. In consideration of the more general nature of DAR over SAR, we skipped SAR and directly applied DAR modeling to reanalyze the hot spring microbiome datasets of Sharp et al. (2014). In fact, our DAR analysis, as presented in previous sections, included SAR as a special case when diversity order q = 0. Our study therefore provides the first glimpse of the SAR/DAR of the hot spring microbiome. The results and conclusions we obtained should certainly be further verified in future with more extensive datasets of the hot spring microbiomes. Although the sample size of 165 hot springs, we used, is not small, the future studies should attempt to collect samples from more diverse regions from different continents to validate our study on a truly global scale.

Author Contributions

ZM designed the study and wrote the paper. LL performed the data analysis. All authors approved the submission.

Conflict of Interest Statement

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

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2019.00118/full#supplementary-material

References

Amin, A., Ahmed, I., Salam, N., Kim, B. Y., Singh, D., Zhi, X. Y., et al. (2017). Diversity and distribution of thermophilic bacteria in hot springs of Pakistan. Microb. Ecol. 74, 116–127. doi: 10.1007/s00248-017-0930-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Arrhenius, O. (1921). Species and area. J. Ecol. 9, 95–99. doi: 10.2307/2255763

CrossRef Full Text | Google Scholar

Barreto, D. P., Conrad, R., Klose, M., Claus, P., and Enrich-Prast, A. (2014). Distance-decay and taxa-area relationships for bacteria, archaea and methanogenic archaea in a tropical lake sediment. PLoS ONE 9:e110128. doi: 10.1371/journal.pone.0110128

PubMed Abstract | CrossRef Full Text | Google Scholar

Bell, T. (2010). Experimental tests of the bacterial distance-decay relationship. ISME J. 4, 1357–1365. doi: 10.1038/ismej.2010.77

PubMed Abstract | CrossRef Full Text | Google Scholar

Bell, T., Ager, D., Song, J., Newman, J. A., Thompson, I. P., Lilley, A. K., et al. (2005). Larger islands house more bacterial taxa. Science 308:1884. doi: 10.1126/science.1111318

PubMed Abstract | CrossRef Full Text | Google Scholar

Caporaso, J. G., Kuczynski, J., Stombaugh, J., Bittinger, K., Bushman, F. D., Costello, E. K., et al. (2010). QIIME allows analysis of high-throughput community sequencing data. Nat. Methods 7, 335–336. doi: 10.1038/nmeth.f.303

PubMed Abstract | CrossRef Full Text | Google Scholar

Chan, C. S., Chan, K. G., Ee, R., Hong, K. W., Urbieta, M. S., Donati, E. R., et al. (2017). Effects of physiochemical factors on prokaryotic biodiversity in Malaysian circumneutral hot springs. Front. Microbiol. 8:1252. doi: 10.3389/fmicb.2017.01252

PubMed Abstract | CrossRef Full Text | Google Scholar

Chao, A., Chiu, C. H., and Hsieh, T. C. (2012). Proposing a resolution to debates on diversity partitioning. Ecology 93, 2037–2051. doi: 10.1890/11-1817.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Chao, A., Chiu, C. H., and Jost, L. (2014). Unifying species diversity, phylogenetic diversity, functional diversity and related similarity and differentiation measures through Hill numbers. Annu. Rev. Ecol. Evol. Syst. 45, 297–324. doi: 10.1146/annurev-ecolsys-120213-091540

CrossRef Full Text | Google Scholar

Cole, J. K., Peacock, J. P., Dodsworth, J. A., Williams, A. J., Thompson, D. B., Dong, H., et al. (2013). Sediment microbial communities in Great Boiling Spring are controlled by temperature and distinct from water communities. ISME J. 7, 718–729. doi: 10.1038/ismej.2012.157

PubMed Abstract | CrossRef Full Text | Google Scholar

De León, K. B., Gerlach, R., Peyton, B. M., and Fields, M. W. (2013). Archaeal and bacterial communities in three alkaline hot springs in Heart Lake Geyser Basin, Yellowstone National Park. Front. Microbiol. 4:330. doi: 10.3389/fmicb.2013.00330

CrossRef Full Text | Google Scholar

de Oliveira, T. B., Gomes, E., and Rodrigues, A. (2015). Thermophilic fungi in the new age of fungal taxonomy. Extremophiles 19, 31–37. doi: 10.1007/s00792-014-0707-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Fierer, N., Leff, J. W., Adams, B. J., Nielsen, U. N., Bates, S. T., Lauber, C. L., et al. (2012). Cross-biome metagenomic analyses of soil microbial communities and their functional attributes. Proc. Natl. Acad. Sci. U.S.A. 109, 21390–21395. doi: 10.1073/pnas.1215210110

PubMed Abstract | CrossRef Full Text | Google Scholar

Gajigan, A. P., Yñiguez, A. T., Villanoy, C. L., San Diego-McGlone, M. L., Jacinto, G. S., and Conaco, C. (2018). Diversity and community structure of marine microbes around the Benham Rise underwater plateau, northeastern Philippines. Peer J. 6:e4781. doi: 10.7717/peerj.4781

PubMed Abstract | CrossRef Full Text | Google Scholar

Geesey, G. G., Barkay, T., and King, S. (2016). Microbes in mercury-enriched geothermal springs in western North America. Sci. Tot. Environ. 569–570, 321–331. doi: 10.1016/j.scitotenv.2016.06.080

CrossRef Full Text | Google Scholar

Ghilamicael, A. M., Budambula, N. L. M., Anami, S. E., Mehari, T., and Boga, H. I. (2017). Evaluation of prokaryotic diversity of five hot springs in Eritrea. BMC Microbiol. 17:203. doi: 10.1186/s12866-017-1113-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Green, J., and Bohannan, B. J. (2006). Spatial scaling of microbial biodiversity. Trends Ecol. Evol. 21, 501–507. doi: 10.1016/j.tree.2006.06.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Green, J. L., Holmes, A. J., Westoby, M., Oliver, I., Briscoe, D., Dangerfield, M., et al. (2004). Spatial scaling of microbial eukaryote diversity. Nature 432, 747–750. doi: 10.1038/nature03034

PubMed Abstract | CrossRef Full Text | Google Scholar

Hartmann, M., Niklaus, P. A., Zimmermann, S., Schmutz, S., Kremer, J., Abarenkov, K., et al. (2014). Resistance and resilience of the forest soil microbiome to logging-associated compaction. ISME J. 8, 226–244. doi: 10.1038/ismej.2013.141

PubMed Abstract | CrossRef Full Text | Google Scholar

Hedlund, B. P., Murugapiran, S. K., Alba, T. W., Levy, A., Dodsworth, J. A., Goertz, G. B., et al. (2015). Uncultivated thermophiles: current status and spotlight on ‘Aigarchaeota’. Curr. Opin. Microbiol. 25, 136–145. doi: 10.1016/j.mib.2015.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Hill, M. O. (1973). Diversity and evenness: a unifying notation and its consequences. Ecology 54, 427–432. doi: 10.2307/1934352

CrossRef Full Text | Google Scholar

Horner-Devin, C., Lage, M. M., Hughes, J. B., and Bohannan, B. J. (2004). A taxa-area relationship for bacteria. Nature 432, 750–753. doi: 10.1038/nature03073

CrossRef Full Text | Google Scholar

Huttenhower, C., Gevers, D., Knight, R., Abubucker, S., Badger, J. H., Chinwalla, A. T., et al. (2012). Structure, function and diversity of the healthy human microbiome. Nature 486, 207–214. doi: 10.1038/nature11234

CrossRef Full Text | Google Scholar

Inskeep, W. P., Jay, Z. J., Herrgard, M. J., Kozubal, M. A., Rusch, D. B., Tringe, S. G., et al. (2013a). Phylogenetic and functional analysis of metagenome sequence from high-temperature archaeal habitats demonstrate Linkages between metabolic potential and geochemistry. Front. Microbiol. 4:95. doi: 10.3389/fmicb.2013.00095

PubMed Abstract | CrossRef Full Text | Google Scholar

Inskeep, W. P., Jay, Z. J., Tringe, S. G., Herrgård, M. J., and Rusch, D. B. (2013b). The YNP metagenome project: environmental parameters responsible for microbial distribution in the Yellowstone geothermal ecosystem. Front. Microbiol. 4:67. doi: 10.3389/fmicb.2013.00067

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, Z., Li, P., Jiang, D., Dai, X., Zhang, R., Wang, Y., et al. (2015). Microbial community structure and arsenic biogeochemistry in an acid vapor-formed spring in Tengchong Geothermal Area, China. PLoS ONE 11:e0146331. doi: 10.1371/journal.pone.0146331

PubMed Abstract | CrossRef Full Text | Google Scholar

Jost, L. (2007). Partitioning diversity into independent alpha and beta components. Ecology 88, 2427–2439. doi: 10.1890/06-1736.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, K. H., Ding, X. W., Salam, N., Zhang, B., Tang, X. F., Deng, B., et al. (2018). Unexpected fungal communities in the Rehai thermal springs of Tengchong influenced by abiotic factors. Extremophiles 22, 525–535. doi: 10.1007/s00792-018-1014-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, Z. S. (2017). Measuring microbiome diversity and similarity with hill numbers. Metagenomics 2018, 157–178 doi: 10.1016/B978-0-08-102268-9.00008-2

CrossRef Full Text | Google Scholar

Ma, Z. S. (2018a). Extending species-area relationships (SAR) to diversity-area relationships (DAR). Ecol. Evol. 8:10023–10038.

Google Scholar

Ma, Z. S. (2018b). Sketching the human microbiome biogeography with DAR (Diversity Area Relationship) profiles. Microb. Ecol.

PubMed Abstract | Google Scholar

Ma, Z. S., and Li, L. W. (2019). Semen microbiome biogeography: an analysis based on a Chinese population study. Front. Microbiol. (accepted). doi: 10.3389/fmicb.2018.03333

CrossRef Full Text | Google Scholar

MacArthur, R. H., and Wilson, E. O. (1967). The Theory of Island Biogeography. Princeton, NJ: Princeton University Press.

Google Scholar

Masaki, Y., Tsutsumi, K., Hirano, S. I., and Okibe, N. (2016). Microbial community profiling of the Chinoike Jigoku (“Blood Pond Hell”) hot spring in Beppu, Japan: isolation and characterization of Fe(III)-reducing Sulfolobus sp. Strain GA1. Res. Microbiol. 167, 595–603. doi: 10.1016/j.resmic.2016.04.011

CrossRef Full Text | Google Scholar

Merkel, A. Y., Pimenov, N. V., Rusanov, I. I., Slobodkin, A. I., Slobodkina, G. B., Tarnovetckii, I. Y., et al. (2017). Microbial diversity and autotrophic activity in Kamchatka hot springs. Extremophiles 21, 307–317. doi: 10.1007/s00792-016-0903-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyer-Dombard, D. R., and Amend, J. P. (2014). Geochemistry and microbial ecology in alkaline hot springs of Ambitle Island, Papua New Guinea. Extremophiles 18, 763–778. doi: 10.1007/s00792-014-0657-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohanrao, M. M., Singh, D. P., Kanika, K., Goyal, E., and Singh, A. K. (2016). Deciphering the microbial diversity of Tattapani water spring using metagenomic approach. Int. J. Agri. Sci. Res. 6, 371–382. doi: 10.1016/j.gdata.2015.04.014

CrossRef Full Text | Google Scholar

Noguez, A. M., Escalante, A. E., García-Oliva, F., Arita, H. T., Forney, L., and Souza, V. (2005). Microbial macroecology: highly structured prokaryotic soil assemblages in a tropical deciduous forest. Glob. Ecol. Biogeogr. 14, 241–248. doi: 10.1111/j.1466-822X.2005.00156.x

CrossRef Full Text | Google Scholar

Oliverio, A. M., Power, J. F., Washburne, A., Cary, S. C., Stott, M. B., and Fierer, N. (2018). The ecology and diversity of microbial eukaryotes in geothermal springs. ISME J. 12, 1918–1928. doi: 10.1038/s41396-018-0104-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Peay, K. G., Bruns, T. D., Kennedy, P. G., Bergemann, S. E., and Garbelotto, M. (2007). A strong species–area relationship for eukaryotic soil microbes: island size matters for ectomycorrhizal fungi. Ecol. Lett. 10, 470–480. doi: 10.1111/j.1461-0248.2007.01035.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Plotkin, J. B., Potts, M. D., Yu, D. W., Bunyavejchewin, S., Condit, R., Foster, R., Hubbell, S., et al. (2000). Predicting species diversity in tropical forests. Proc. Natl. Acad. Sci. U.S.A. 97, 10850–10854. doi: 10.1073/pnas.97.20.10850

PubMed Abstract | CrossRef Full Text | Google Scholar

Poddar, A., and Das, S. K. (2018). Microbiological studies of hot springs in India: a review. Arch. Microbiol. 200, 1–18. doi: 10.1007/s00203-017-1429-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Pop Ristova, P. P., Wenzhöfer, F., Ramette, A., Felden, J., and Boetius, A. (2014). Spatial scales of bacterial community diversity at cold seeps (Eastern Mediterranean Sea). ISME J. 9, 1306–1318. doi: 10.1038/ismej.2014.217

PubMed Abstract | CrossRef Full Text | Google Scholar

Preston, F. W. (1960). Time and space and the variation of species. Ecology 41, 611–627. doi: 10.2307/1931793

CrossRef Full Text | Google Scholar

Preston, F. W. (1962). The canonical distribution of commonness and rarity: Part, I. Ecology 43, 185–215. doi: 10.2307/1931976

CrossRef Full Text | Google Scholar

Ruff, S. E., Biddle, J. F., Teske, A. P., Knittel, K., Boetius, A., and Ramette, A. et al. (2015). Global dispersion and local diversification of the methane seep microbiome. Proc. Natl. Acad. Sci. U.S.A. 112, 4015–4020. doi: 10.1073/pnas.1421865112

PubMed Abstract | CrossRef Full Text | Google Scholar

Salano, O. A., Makonde, H. M., Kasili, R. W., Wangai, L. N., Nawiri, M., and Boga, H. (2017). Diversity and distribution of fungal communities within the hot springs of soda lakes in the Kenyan rift valley. Afr. J. Microbiol. Res. 19, 764–777. doi: 10.5897/AJMR2017.8510

CrossRef Full Text | Google Scholar

Schloss, P. D., Westcott, S. L., Ryabin, T., Hall, J. R., Hartmann, M., Hollister, E. B., et al. (2009). Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl. Environ. Microbiol. 75:7537. doi: 10.1128/AEM.01541-09

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharp, C. E., Brady, A. L., Sharp, G. H., Grasby, S. E., Stott, M. B., and Dunfield, P. F. (2014). Humboldt's spa: microbial diversity is controlled by temperature in geothermal environments. ISME J. 8, 1166–1174. doi: 10.1038/ismej.2013.237

PubMed Abstract | CrossRef Full Text | Google Scholar

Sinha, R., Abnet, C. C., White, O., Knight, R., and Huttenhower, C. (2015). The microbiome quality control project: baseline study design and future directions. Genome Biol. 16:276. doi: 10.1186/s13059-015-0841-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, V. H., Foster, B. L., Grover, J. P., Holt, R. D., Leibold, M.A., and deNoyelles, F. (2005). Phytoplankton species richness scales consistently from laboratory microcosms to the world's oceans. Proc. Natl. Acad. Sci. U.S.A. 102, 4393–4396. doi: 10.1073/pnas.0500094102

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, Z. Q., Wang, F. P., Zhi, X. Y., Chen, J. Q., Zhou, E. M., Liang, F., et al. (2013). Bacterial and archaeal diversities in Yunnan and Tibetan hot springs, China. Environ. Microbiol. 15, 1160–1175. doi: 10.1111/1462-2920.12025

PubMed Abstract | CrossRef Full Text | Google Scholar

Terrat, S., Dequiedt, S., Horrigue, W., Lelievre, M., Cruaud, C., Saby, N. P. A., et al. (2015). Improving soil bacterial taxa–area relationships assessment using DNA meta-barcoding. Heredity 114, 468–475. doi: 10.1038/hdy.2014.91

PubMed Abstract | CrossRef Full Text | Google Scholar

Tjørve, E. (2009). Shapes and functions of species–area curves (II): a review of new models and parameterizations. J. Biogeogr. 36, 1435–1445. doi: 10.1111/j.1365-2699.2009.02101.x

CrossRef Full Text | Google Scholar

Ulrich, W., and Buszko, J. (2003). Self-similarity and the species–area relation of Polish butterflies. Basic Appl. Ecol. 4, 263–270. doi: 10.1078/1439-1791-00139

CrossRef Full Text | Google Scholar

Urbieta, M. S., Donati, E. R., Chan, K., Shahar, S., Sin, L. L., and Goh, K. M. (2014a). Thermophiles in the genomic era: biodiversity, science, and applications. Biotechnol. Adv. 33, 633–647. doi: 10.1016/j.biotechadv.2015.04.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Urbieta, M. S., Toril, E. G., Giaveno, M. A., Bazán, A. A., and Donati, E. R. (2014b). Archaea land bacterial diversity in five different hydrothermal ponds in the Copahue region in Argentina. Syst. Appl. Microbiol. 37, 429–441. doi: 10.1016/j.syapm.2014.05.012

CrossRef Full Text | Google Scholar

Várbíró, G., Görgényi, J., Tóthmérész, B., Padisák, J., Hajnal, E., and Borics, G. (2017). Functional redundancy modifies species–area relationship for freshwater phytoplankton. Ecol. Evol. 7, 9905–9913. doi: 10.1002/ece3.3512

PubMed Abstract | CrossRef Full Text | Google Scholar

Watson, H. C. (1835). Remarks on Geographic Distribution of British Plants. London: Longman; Rees; Orme; Brown; Green and Longman.

Xie, W., Zhang, C. L., Wang, J., Chen, Y., Zhu, Y., de la Torre, J. R., et al. (2015). Distribution of ether lipids and composition of the archaeal community in terrestrial geothermal springs: impact of environmental variables. Environ. Microbiol. 17, 1600–1614. doi: 10.1111/1462-2920.12595

PubMed Abstract | CrossRef Full Text | Google Scholar