Bacterial diversity in a glacier foreland of the high Arctic
Ursel Schütte is interested in understanding patterns of biogeography for microorganisms by applying macroecological concepts to the microbial world. Zaid Abdo is a professor of Statistics and Mathematics and he combines mathematical modeling and statistical analysis using computationally intensive methods and simulations to address current problems in biology such as the analysis of microbial community composition and diversity, experimental evolution, and the barcoding of life. James Foster explores and attempts to understand both natural and simulated evolution by developing and analyzing genetic algorithms. He also investigates biological problems such as the emergence and structure of microbial ecosystems. Jacques Ravel's research interests extend from the application of microbial genomics to explore the human microbiome, to comparative analyses of microbial genome sequences and microbial forensics. John Bunge's research focuses on the problem of estimating the species richness or biodiversity in a community. Bjørn Solheim focuses on Arctic terrestrial ecosystems, in particular cyanobacterial communities and their role in nitrogen fixation, and the impact of climate change on these types of communities. Larry Forney's research centers on the diversity and distribution of prokaryotes. Both field and laboratory studies are done to explore the temporal and spatial patterns of community diversity, as well as factors that influence the dynamics of inter- and intra-species competition and how environmental conditions might influence the tempo of adaptive evolution.
Abstract
Over the past 100 years, Arctic temperatures have increased at almost twice the global average rate. One consequence is the acceleration of glacier retreat, exposing new habitats that are colonized by microorganisms whose diversity and function are unknown. Here, we characterized bacterial diversity along two approximately parallel chronosequences in an Arctic glacier forefield that span six time points following glacier retreat. We assessed changes in phylotype richness, evenness and turnover rate through the analysis of 16S rRNA gene sequences recovered from 52 samples taken from surface layers along the chronosequences. An average of 4500 sequences was obtained from each sample by 454 pyrosequencing. Using parametric methods, it was estimated that bacterial phylotype richness was high, and that it increased significantly from an average of 4000 (at a threshold of 97% sequence similarity) at locations exposed for 5 years to an average of 7050 phylotypes per 0.5 g of soil at sites that had been exposed for 150 years. Phylotype evenness also increased over time, with an evenness of 0.74 for 150 years since glacier retreat reflecting large proportions of rare phylotypes. The bacterial species turnover rate was especially high between sites exposed for 5 and 19 years. The level of bacterial diversity present in this High Arctic glacier foreland was comparable with that found in temperate and tropical soils, raising the question whether global patterns of bacterial species diversity parallel that of plants and animals, which have been found to form a latitudinal gradient and be lower in polar regions compared with the tropics.
Introduction
The Arctic is disproportionately affected by climate change and arctic temperatures have increased at almost twice the global average rate over the past 100 years (Bernstein et al. 2007). The acceleration of glacier retreat is one consequence of these temperature increases and the net mass balance of most Arctic glaciers has been negative over the past few decades (Dowdeswell et al. 1997). Predictions are that the rate of loss will increase in the future (Kintisch 2009). The newly exposed rocks and minerals represent new habitats for microorganisms that have key roles in soil development, nutrient cycling and plant growth. Bacterial communities may be key determinants of Arctic ecosystem stability and function because of their important roles in soil development, biogeochemical cycles and heterotrophic activities. As fundamental knowledge on bacterial communities in High Arctic glacier forelands is lacking, studies on bacterial community diversity and composition in these extreme environments are needed. The design of these studies should take into account the extraordinary spatial heterogeneity of these habitats and patchy distribution of bacteria in soils (Green et al. 2004; Horner-Devine et al. 2004; Green & Bohannan 2006).
Valuable information has been gained from studies carried out to assess bacterial diversity and succession in montane glacier forelands (Schipper et al. 2001, Sigler & Zeyer 2002; Sigler et al. 2002, Nemergut et al. 2007; Schmidt Reed et al. 2008). These studies that were based on the analysis of 16S rRNA gene sequences provide conflicting information pertaining to changes in phylotype diversity over time following glacier retreat. Nemergut et al. (2007) found that it increased while others found the opposite (Sigler & Zeyer 2002; Sigler et al. 2002). These findings suggest that the development of bacterial communities in soils following glacier retreat may be less predictable than in plant communities, where species richness consistently increases over time (Huston & Smith 1987; Matthews 1992; Hodkinson et al. 2003; Walker & del Moral 2003). Although the methods used were appropriate for the aims of these studies, the sampling schemes, numbers of samples and the use of low resolution analytical methods (Bent & Forney 2008) precluded statistically meaningful estimates of bacterial species richness (α-diversity), evenness or turnover rate (β-diversity). Furthermore, the extrapolation of findings from studies of montane glacier forelands to those of the High Arctic may be problematic because the dynamics of alpine and High Arctic glaciers differ. High Arctic glaciers are commonly polythermal and they gather a large reservoir of meltwater beneath the glacial mass during winter, which is then released abruptly as a surging breakout in spring. These thermal and hydrological differences result in comparatively unstable forelands in the High Arctic (Harland 1997; Hodkinson et al. 2003). Thus, fundamental differences might exist between the montane and High Arctic glacier forelands.
The few studies carried out provide only limited insight to microbial diversity in polar soils as only a small number of samples were examined. Nonetheless, the findings suggest that bacterial diversity in polar habitats is high despite the extreme environmental conditions that exist (Broady 1996; Vishniac 1996, Shi et al. 1997; Zhou et al. 1997; Lawley et al. 2004, Shivaji et al. 2004; Neufeld & Mohn 2005; Steven et al. 2007; González-Toril et al. 2008), and a large proportion of the bacterial taxa are unlike those previously described. An analysis of Siberian tundra soil showed that mainly α-, δ-Proteobacteria were present, but 7% of the 16S rRNA gene sequences retrieved from the samples were more than 20% different from those previously reported (Zhou et al. 1997). The findings were in agreement with those of Shivaji et al. (2004) who found members of seven different phyla including α-, β-, γ-Proteobacteria and Actinobacteria, and culturable bacteria that belonged to genera previously not reported to be in Antarctic soils such as Sphingomonas and Acidovorax. Neufeld & Mohn (2005) did probably the most extensive study on bacterial diversity in this environment. The authors sampled geographically distant boreal and arctic soils, and diversity was studied using ∼13 000 serial analysis of ribosomal sequence tags and denaturing gradient gel electrophoresis of 16S rRNA genes. They reported that diversity in arctic soils was higher than in boreal soils with the highest diversity found at 82 °N, an extreme northern location. Until recently, investigations of bacterial diversity based on cultivation-independent methods were limited by the lack of high throughput methods for sequencing 16S rRNA genes. The advent of ‘next generation’ DNA sequencing technologies such as pyrosequencing has obviated this problem and provided the means necessary to conduct intensive and extensive analyses of bacterial diversity in different environments (Sogin et al. 2006; Roesch et al. 2007; Dethlefsen et al. 2008; Huse et al. 2008).
We recently showed that bacterial succession occurred in the glacier foreland of Midtre Loven glacier (Schütte et al. 2009), and here we report the results of studies carried out to assess species (phylotype) richness, evenness and turnover rate along two chronosequences in this foreland. These approximately parallel chronosequences span six times of exposure following retreat of the Midtre Lovén glacier. The bacterial diversity in surface soils was assessed by pyrosequencing of bar-coded 16S rRNA gene amplicons present in samples taken from surface soils. We did not expect high diversity in glacier foreland soils because environmental conditions are harsh including UV radiation, low water content and low nutrient contents. Furthermore, increase in plant richness and higher soil organic matter and higher water content (Hodkinson et al. 2003) most probably results in an increase in the number of spatial niches in the soil thus, leading to the expectation that bacterial phylotype richness increases over time in the foreland of Midtre Lovén glacier. Using different statistical methods to estimate species richness, evenness and the species turnover rate, we were able to show a remarkably high bacterial diversity at all times since glacier retreat with an increase over time.
Materials and methods
Study site and sampling
The foreland of Midre Lovén glacier near the settlement Ny-Ålesund, West Spitsbergen (74–81 °N; 10–35 °E) was chosen as a field site to assess bacterial diversity. Samples were taken in the beginning of August 2003 along the chronosequence established by Hodkinson et al. (2003) and along a second approximately parallel chronosequence located about 25 m away. The second chronosequence was designed to account for differences in bacterial community structure caused by horizontal variation and the two chronosequences were assumed to be independent of each other (Schütte et al. 2009). Both chronosequences included six sites with a total distance of about 1 km and the distance between sites was approximately 100–200 m. The total time period covered was 150 years with site 1–6 representing times since glacier retreat of 5, 19, 40, 63, 100 and 150 years respectively (Schütte et al. 2009). Hodkinson et al. (2003) had previously determined the time since deglaciation of each sampling location using photographs and radiocarbon analysis.
At each sampling location, a transect 5 m long was established along which five samples were taken 1 m apart. Each sample had a surface area of 10 × 10 cm2. Larger gravel (> ∼3–5 cm) was manually removed before the samples were taken. In this study, only the samples taken from the surface were used, which contained a mixture of vegetation, rhizosphere and bulk soil and did not vary in thickness among sampling sites (about 1 cm). In total, 60 samples were collected, placed in plastic bags and kept frozen on ice during transport to the laboratory and stored at −80 °C until analysed.
454 Pyrosequencing
Genomic DNA was isolated from 0.5 g of the soil samples, which had been mixed in their bag, using a modified procedure based on the UltraClean™ Fecal DNA Kit (MoBio). For each sample, soil was weighed into a sterile 2 mL tube containing 0.5 g of silicate glass beads (Glenn Mills) and 750 μL of TE (1 mM Tris and 50 mM EDTA) was added. The cells were treated with 50 μL lysozyme (10 mg/mL; Sigma-Aldrich) and 25 μL of mutanolysin (2 mg/mL; Sigma-Aldrich). Afterwards, 60 μL of the solution S1 and 200 μL of the IRS solution (MoBio) were added followed by a bead beating step for 3 min at full speed (Biospec Products Inc). The tubes were centrifuged for 4 min at 13 000 g and 450 μL of the supernatants was transferred to a new 2 mL tube before following the protocol from MoBio. The columns were incubated for 5 min at room temperature before eluting the DNA using 50 μL deionized water with pH 7, heated to 90 °C. As positive controls, we used cells from pure cultures of Escherichia coli K12 MG1655 and Lactobacillus aviaris ATCC 43232. The DNA from all samples was further cleaned using the Wizard PCR purification kit (Promega).
We amplified the bacterial 16S rRNA gene and the primers used flanked the variable regions 1 and 2 (Escherichia coli positions 27F-338R). The primer sequences were as follows: 454_27F 5′-GCCTTGCCAGCCCGCTCAGTCAGAGTTTGATCCTGGCTCAG and 454_338R 5′-GCCTCCCTCGCGCCATCAGNNNNNNNNCATGCTGCCTCCCGTAGGAGT; where the underlined sequences are 454 Life Sciences® adapters B and A in 27F and 338R respectively, and the bold font denotes the universal 16S rRNA primers 27F and 338R. The 338R primer included a unique sequence tag to barcode each of the samples denoted by the eight N. This allowed us to sequence the amplicons from all samples simultaneously, and afterwards assign each sequence to the sample they were obtained from. Using a specific barcode for each reverse primer resulted in each PCR reaction containing a different reverse primer, thus, a negative control was run for each reaction. Each PCR reaction contained 35.2 μL of PCR-grade water, 5.0 μL of 10× buffer (Applied Biosystems), 6.0 μL of 25 mM MgCl2 (Applied Biosystems), 1.2 μL of 20 mg/mL BSA (Roche Applied Science), 0.4 μL of 25 mM dNTP (Amersham Bioscience), 0.5 μL of 20 μM forward primer 454_27f, 1 μL of 10 μM reverse primer 454_338r, 0.2 μL of 5 U/μL Taq DNA polymerase (Applied Biosystems) and 1.0 μL of DNA template, in a total volume of 50 μL. Amplification of 16S rRNA genes was performed using an initial denaturation step at 94 °C for 4 min, followed by 34 cycles of denaturation at 94 °C for 1 min, annealing at 55 °C for 1 min and an extension at 72 °C for 2 min. The final extension was 10 min step at 72 °C. The concentrations of amplicons were estimated using a GelDoc quantification system (Bio-Rad), and approximately equal amounts (∼100 ng) of all amplicons were mixed in a single tube. Amplification primers and reaction buffer were removed by processing the amplicon mixture with the AMPure Kit (Agencourt). Emulsion PCRs were performed as described in Margulies et al. (2005). Sequences were obtained using a Roche 454 GS-FLX (Roche-454 Life Sciences).
Quality control of sequencing
We followed the criteria previously described to assess the quality of sequence reads (Sogin et al. 2006; McKenna et al. 2008). To pass, a sequence read had to (i) include a perfect match to the sequence tag (bar-code) and the 16S rRNA gene primer; (ii) be at least 100 bp in length; (iii) have no more than four undetermined bases; and (iv) have at least 60% match to a previously determined 16S rRNA gene sequence after alignment with NAST (DeSantis et al. 2006). Sequence reads with the same bar-code were grouped using a script that integrates the whole genome alignment tool MUMmer (Kurtz et al. 2004) with the 16S rRNA gene reverse primer and bar-code sequences as references.
The error rate in pyrosequencing is reportedly 0.25% per read (Huse et al. 2007), suggesting that most of the sequence reads obtained after this filtering contained either 0 or 1 error. As demonstrated previously (McKenna et al. 2008), this potential single nucleotide error will not interfere with data analysis and is unlikely to have caused a sequence to be assigned to the incorrect taxon because more than 12–14 errors would be required to form a new taxon at the species level using an operational taxonomic unit threshold of 97% sequence similarity (McKenna et al. 2008). The assignment of sequences to taxonomic units was carried out using the Ribosomal Database Project Bayesian classifier (Wang et al. 2007).
Statistical analysis
Estimation of species richness. To investigate changes in phylotype richness along the chronosequences, both parametric and nonparametric richness estimation methods were used. This was carried out for all 52 individual samples taken along the chronosequence. Sequences were aligned within each sample using the secondary-structure alignment program INFERNAL that is accessible through the RDP pyrosequencing webpage (Cole et al. 2009). We then computed evolutionary distances between all sequences within each sample using DNAdist from the Phylip 3.68 package (Felsenstein 1993), with default parameters (F84 distance correction, transition/transversion ratio of 2.0). Given these distances, we used DOTUR (Schloss & Handelsman 2005) to cluster all sequences using furthest-neighbour clustering with a 1%, 3% and 10% similarity thresholds. The resulting abundance distributions were then converted into frequency counts to determine how many sequences occurred once (singletons), twice (doubletons) and so on. The approaches for parametric richness estimation were previously described by Bunge & Barger (2008). In brief, different models were fitted to the frequency counts for all possible truncation points of the frequency counts. Different goodness of fit tests (GF0–GF5) and the Akaike Information Criterion (AIC) were then used to determine how well each model fit the frequency count data (see Bunge & Barger (2008) for details). For a model to be chosen to fit the data best, GF0–GF5 had to be greater than 0.01, and the AIC had to be lower than that for other models. We used estimates based on finite mixture models, because these models performed generally better than others (see Bunge & Barger (2008) for details). For example, the Inverse Gaussian model frequently resulted in unreasonably large estimates with large standard errors. For nonparametric estimators, the ACE and ACE1 estimates of richness were used (Chao & Bunge 2002). If the coefficient of variation (CV) of ACE was greater than 0.8, the estimate from ACE1 was used and if CV less than 0.8, the estimate from ACE was used (Chao & Bunge 2002). Curves were fitted to the parametric and nonparametric richness estimates observed for each different threshold of sequence similarity to describe the change in richness over time. The significance of the trend was tested using analysis of variance (anova) based on stepwise model selection starting with the most parameter-rich model (fourth order model). Significance of the model was then determined correcting for multiple testing. To ensure that richness estimation was not influenced by low number of sequences per individual sample, the analysis was repeated, and this time the sequences from all the replicates for the corresponding time of exposure since glacier retreat were pooled.
Calculation of evenness. The evenness index Evar (Smith & Wilson 1996) was used to quantify evenness at different locations along the chronosequences. Evar is mathematically independent from richness and it has the property of being symmetric, that is, it weighs abundant and rare taxa equally (Beisel et al. 2003). Although this index may have an overall lower sensitivity (Beisel et al. 2003), symmetry is an essential characteristic, particularly when measuring evenness in bacterial communities that have large numbers of rare members.
Calculation of turnover rate. The changes in species composition over time (species turnover rate or beta diversity) along the chronosequences were calculated using a two-step procedure because it was computationally challenging to determine pairwise distances for all 230 000 sequences. First, all sequences from the time of exposure were pooled and pairwise distances between sequences were calculated using DNAdist (F84 correction with default transition transversion ratio, Felsenstein 1993). Sequences were then clustered based on a threshold of 97% sequence similarity using the furthest-neighbour algorithm (DOTUR, Schloss & Handelsman 2005). Second, for each cluster, we selected the first FASTA sequence identifier as a proxy. These sequences representing the different clusters from the six times since glacier retreat were then pooled. To determine the similarity of these sequences, all pairwise distances were calculated using DNAdist (Felsenstein 1993), and the sequences were clustered into phylotypes based on a threshold of 97% sequence similarity using DOTUR (Schloss & Handelsman 2005). After the clustering, a table was made that showed which phylotype was present at which time since glacier retreat. Based on this matrix, the turnover rate was then calculated using the measures described by Whittaker (1972) and Diamond & May (1977) using the R package (R Foundation for Statistical Computing).
Results
In this study, we assessed bacterial diversity along two approximately parallel chronosequences in a glacier foreland in the High Arctic based on phylotype richness, evenness and turnover rate. This is challenging because different methods use different metrics that lead to disparities in the estimates. For example, nonparametric methods to estimate richness consistently provide lower estimates of richness than do parametric methods (Bunge & Barger 2008). Here, we used multiple approaches for each aspect of diversity investigated and these provided a range of values for phylotype richness, evenness and turnover rates. Various thresholds of sequence similarity among 16S rRNA gene sequences are commonly used as a proxy for different taxonomic levels in studies on microbial diversity with 90%, 97% and 99% corresponding to differences between genera, species and strains respectively (Hagström et al. 2002, Oakley et al. 2008).
Changes in phylotype richness
Phylotype richness was estimated and changes were determined along the chronosequences using both nonparametric (ACE and ACE1, Chao & Lee 1992) and parametric richness estimates (mixture models, Bunge & Barger 2008). Estimates of phylotype richness were obtained for three thresholds of sequence similarity for the replicate soil samples taken from each time since glacier retreat (Fig. 1). Richness estimates obtained from the corresponding samples from two chronosequences were in agreement (data not shown). For all times since glacier retreat, nonparametric estimates were consistently lower than parametric estimates, but both approaches indicated that phylotype richness was high. For example, using a threshold of 97% sequence similarity phylotype richness at the different locations ranged from 2700 to 5800 per 0.5 g of soil, whereas the parametric estimates ranged on average from 4100 to 7000 per 0.5 g of soil. As expected, estimates of richness changed as the threshold for sequence similarity varied. For example, at a threshold of 90% sequence similarity, there was an average of 2700 phylotypes per 0.5 g of soil at 5 years since glacier retreat (nonparametric estimate), but with a threshold of 99% sequence similarity, the richness estimate averaged 6000 phylotypes per 0.5 g of soil (nonparametric estimation, Fig. 1). Furthermore, phylotype richness varied among replicates within each time of glacier retreat (Fig. 1). For example, at 63 years after glacier retreat, the lowest estimate for richness at 97% sequence similarity was approximately 4000 phylotypes per 0.5 g of soil and the largest was about 7400 phylotypes per 0.5 g of soil (parametric estimation). Using regression analysis in conjunction with anova based on stepwise model selection (data not shown), we found that phylotype richness increased significantly over time despite spatial heterogeneity (Fig. 1). This was the case for all three thresholds of sequence similarity, which suggests that phylotype richness changed not only on a fine taxonomic scale (99% sequence similarity), but also on a broader taxonomic scale (90% sequence similarity). As the number of sequences in each sample was low, the standard errors, not among the replicates, but obtained when fitting the parametric and nonparametric models to the frequency count data, were large for some of the individual richness estimates. By pooling the sequences from all the replicates of the corresponding time since glacier retreat, we had about 35 000 sequences per sample, standard errors for each of the richness estimates were low and we were able to confirm that phylotype richness increased over time (data not shown). Using this approach, we were able to show that for all three thresholds of sequence similarity, phylotype richness was high and increased significantly over time independent of the method used to estimate richness.
Estimation of phylotype richness of bacterial communities sampled from six times since glacier retreat along the chronosequence based on data from 454 pyrosequencing of 16S rRNA genes. Richness estimation was carried out for three levels of sequence similarities (90%, 97% and 99%) using nonparametric and parametric richness estimators (Bunge & Barger 2008). Estimates from the corresponding samples from the two chronosequences were in good agreement and thus were not shown separately. If the coefficient of variance (CV) of the frequency data was greater than 0.8, the estimate ACE1 was used; otherwise ACE was used. For a parametric model (either two-mixed exponential or three-mixed exponential) to be chosen to fit the data best, GF0–GF5 had to be greater than 0.01 and the AIC had to be the lowest compared with the other models for that particular truncation point. Regression analysis was performed for each percentage sequence similarity to describe the shape of change in richness over time. The significance of the trend was tested using analysis of variance (anova) based on stepwise model selection starting with the most parameter-rich model and correcting for multiple testing. Linear models were significant except for nonparametric estimation based on 90% sequence similarity in which case a third order model was significant.
Changes in phylotype evenness
The structure of a bacterial community is characterized also by the evenness of the phylotype rank-abundance distribution. Thus, we determined changes in evenness along the chronosequence in the foreland of Midre Lovén glacier using the index Evar. This index is mathematically independent from richness (Smith & Wilson 1996) and weighs variation in rare phylotypes equally to variation in abundant members. This characteristic of symmetry seemed to be essential for determining evenness of these communities given that no phylotypes in these bacterial communities were abundant (about 9% was ‘abundant’ in our samples, but 40% is considered abundant in Smith & Wilson (1996) and Beisel et al. (2003)). As for phylotype richness, the evenness values were determined for sequence similarity thresholds of 90%, 97% and 99% for each of the replicate soil samples taken from each location along the chronosequences (Fig. 2). Evenness did not differ between the two chronosequences (data not shown). Overall, evenness was high considering that log-normal distributions of bacterial species are thought to be common for bacterial communities in natural environments (Curtis et al. 2002). For example, at 150 years since glacier retreat, the evenness values of Evar for 99% sequence similarity were on average 0.74. Evenness increased with finer taxonomic resolution. This is illustrated by the values obtained at 5 years since glacier retreat; the evenness values of Evar for 90% sequence similarity were on average 0.39, at 97% sequence similarity the average was 0.55 and evenness increased to an average of 0.64 at 99% sequence similarity (Fig. 2). This suggests that phylotypes that were more dominant than other phylotypes on a broader taxonomic scale consisted of a number of relatively evenly distributed ‘species’ or ‘strains’, which resulted in an overall higher evenness at a finer taxonomic scale. As for phylotype richness, there was some within-site variation in evenness for each time since glacier retreat. For example, at 40 years since glacier retreat and a threshold of 90% sequence similarity, evenness values of Evar of the individual replicates ranged from 0.4 to above 0.5 (Fig. 2). Evenness values were plotted over time (Fig. 2), but curves were not fitted to the data because the assumption of homoscedasticity (also known as constant variance) was violated even after data transformation. However, although the absolute values of the two evenness indices differ, they are consistent in that both suggest an increase in evenness over time.
Phylotype evenness of bacterial communities sampled from six times since glacier retreat along the chronosequence based on data from 454 pyrosequencing of 16S rRNA genes. Evenness was determined for three levels of sequence similarities (90%, 97% and 99%) using the E1/D and Evar index respectively. Estimates from the corresponding samples from the two chronosequences were in good agreement and thus were not shown separately.
The trends in 1, 2 suggest a positive correlation between phylotype richness and evenness. To assess whether phylotype evenness was significantly correlated with richness, Pearson correlation analysis was carried out for both combinations of richness estimators and evenness measure used in this study. Phylotype richness was significantly positively correlated with evenness for all comparisons and for all thresholds of sequence similarities (Table 1). This positive correlation suggests that in these soils, phylotype richness and evenness may both be meaningful descriptors of bacterial diversity.
| Sequence similarity (%) | Type of richness estimation | Evenness |
|---|---|---|
| E var | ||
| 90 | Nonparametric | 0.624* |
| Parametric | 0.528* | |
| 97 | Nonparametric | 0.718* |
| Parametric | 0.537* | |
| 99 | Nonparametric | 0.700* |
| Parametric | 0.622* |
- *P < 0.001.
Phylotype turnover rates
To measure changes in compositions of bacterial communities along the chronosequence, we determined the turnover rate using the measure by Whittaker (1972) and Diamond & May (1977) based on a 97% sequence similarity. These indices are measures of turnover and they describe the dissimilarity of community compositions. The measure by Diamond & May (1977) additionally accounts for time and thus, it represents a rate of change. Possible values range from zero to one and higher values suggest larger turnover. As a result of computational complexity, these calculations were only carried out for the data sets in which sequences of all replicate soil samples from each time since glacier retreat were pooled. Based on the assumption of the ‘chronosequence approach’ that space can be substituted for time (Walker & del Moral 2003), the results based on the measure by Whittaker (1972) suggest a high turnover in community composition for all pairwise comparisons of bacterial communities along the chronosequence (Table 3). For example, the turnover in phylotypes between 5 and 19 years since glacier retreat suggests that about 98% of the phylotypes between the communities differed. This value was similar to the comparison of the community composition obtained from 5 and 150 years since glacier retreat (0.96, Table 2). The analysis of turnover using the Whittaker index did not take into account that the difference in times since glacier retreat was not constant among all sites sampled. For example, the difference in time since glacier retreat between the first and second location was 14 years, but between locations five and six it was 50 years, and from the first location to the location exposed the longest it was 145 years. Using the method by Diamond & May (1977), the results showed that in the first 14 years turnover rate was about 7% per year (0.071, Table 3) and that this rate decreased to about 0.7% per year from 100 to 150 years since glacier retreat (0.006, Table 3). Overall, these results suggest that there was tremendous turnover of bacterial communities along the chronosequence and that changes in bacterial community compositions were particularly rapid in the early times following glacier retreat.
| Years since glacier retreat | 5 | 19 | 40 | 63 | 100 | 150 |
|---|---|---|---|---|---|---|
| 5 | 0 | 0.071 | 0.029 | 0.017 | 0.010 | 0.007 |
| 19 | 0 | 0.047 | 0.023 | 0.012 | 0.007 | |
| 40 | 0 | 0.043 | 0.016 | 0.009 | ||
| 63 | 0 | 0.027 | 0.011 | |||
| 100 | 0 | 0.02 | ||||
| 150 | 0 |
| Years since glacier retreat | 5 | 19 | 40 | 63 | 100 | 150 |
|---|---|---|---|---|---|---|
| 5 | 0 | 0.978 | 0.978 | 0.97 | 0.984 | 0.96 |
| 19 | 0 | 0.986 | 0.982 | 0.988 | 0.976 | |
| 40 | 0 | 0.977 | 0.977 | 0.973 | ||
| 63 | 0 | 0.981 | 0.964 | |||
| 100 | 0 | 0.977 | ||||
| 150 | 0 |
Bacterial community composition
Similar to the study by Roesch et al. (2007), a large proportion of the 16S rRNA gene sequences was not similar to anything previously reported; 15% of the sequences could not be assigned to any known phylum and on average 70% of the sequences in a sample could not be classified into a genus. The sequences that could be classified were assigned to 20 phyla (Acidobacteria, Chlamydiae, BRC1, Nitrospira, OD1, Chloroflexi, WS3, Deinococcus-Thermus, Bacteroidetes, Proteobacteria, OP11, Firmicutes, Planctomycetes, Spirochaetes, OP10, TM7, Gemmatimonadetes, Actinobacteria, Cyanobacteria, Verrucomicrobia) and 284 genera (see Tables S1 and S2, Supporting information). The 10 most abundant genera were GP4, Spirosoma, Sphingomonas, GP6, GP3, Terromonas, Hymenobacter, Gemmatimonas, Brevundimonas and Sphingopyxis. Many sequences found belonged to genera that grow oligotrophically and chemoorganotrophically such as 19 genera belonging to the phylum Acidobacteria and 54 genera belonging to the Actinobacteria. Diverse genera of nitrogen fixing genera were present including the cyanobacterial genera GPVI, GPIIa, GPI and GPIV, as well as Rhizobium, Mesorhizobium, Bradorhizobium, Hyrdogenophaga, Burkholderia, Herbaspirillum, Azohydromonas, Blastobacter, Methylocystis, Agromonas, Methylosinus, Ensifer, Acetobacterium, Pelomonas, Devosia and Methylocella. The large proportion of the unassigned sequences underlines the large diversity and how little we know about the bacterial diversity in these High Arctic soils.
Discussion
This is the first time that bacterial diversity in a terrestrial ecosystem in the High Arctic has been intensively and extensively assessed using cultivation-independent methods. Replicate soil samples from six locations along two parallel chronosequences were collected in the foreland of Midtre Lovén glacier that span a 150 years period since deglaciation. The diversity was characterized by phylogenetic analysis of 16S rRNA sequences that had been amplified from genomic DNA extracted from the soil samples. By using high throughput pyrosequencing methods, about 35 000 sequences were obtained for each of the six times following deglaciation. Despite the spatial heterogeneity of these habitats, the experimental design and statistical methods used allowed us to demonstrate that bacterial phylotype richness and evenness were higher than we expected even at sites that had been exposed for as little as 5 years and increased more or less linearly over subsequent years.
There was a statistically significant positive correlation between changes in phylotype richness and evenness along the chronosequence. This positive correlation is consistent with the models describing relationships among various diversity indices by deBenedictis (1973) and Hill (1973), but stand in contrast to multiple studies that found either no correlation or a negative relationship between species richness and evenness (Stirling & Wisley 2001; Ma 2005; Bock et al. 2007). The index used in this study to calculate evenness is, in contrast to other indices such as NHC, mathematically independent from phylotype richness (Smith & Wilson 1996). The bacterial communities present in soils of the Midtre Lovén glacier foreland showed high evenness with no dominant members and a large proportion of rare taxa. Using a threshold of 97% sequence similarity, about 25% of the sequences in each sample occurred only once or twice in a sample, while common types such as Acidobacteria only accounted for about 9% of all the sequences. Thus, the increase in species richness reflects the addition of uncommon phylotypes, resulting in a more flat rank-abundance distribution of phylotypes (higher evenness) and thus, a positive correlation between richness and evenness. However, a sampling bias may also partly account for this correlation (Bent & Forney 2008). To illustrate this, consider two bacterial communities that have comparable phylotype richness, but in one the phylotypes are equally abundant, while in the other the rank abundance of phylotypes is skewed with one taxon accounting for 90% of all the members, and additional rare populations accounting for the remaining 10%. DNA sequencing constitutes a ‘random draw’ from the pool of 16S rRNA gene sequences in a sample (Bent & Forney 2008). Hence, in the latter community it is more probably that sequences from the dominant population will be sampled, while in the former community there is a higher probability that a ‘draw’ will be from a taxon not previously sampled and the apparent richness of the community would be higher (Bent & Forney 2008). Fortunately, the sampling bias in this study may be of less importance as the distribution of all populations in these bacterial communities was relatively even.
Perhaps the most striking finding from this study was the extraordinarily high species turnover (beta diversity) that occurred soon after glacier retreat and that continued thereafter. This was consistent with our previous findings that succession occurred along the chronosequences (Schütte et al. 2009). Although this result is based on the assumption of the ‘chronosequence approach’ that space can be substituted for time (Walker & del Moral 2003) and the cause of this species turnover is of course unknown, it is tempting to speculate that it was driven by substantial changes in the soil environment following deglaciation. Alternatively, it could reflect the fact that bacterial populations initially present were not well suited to the extant conditions in newly exposed soils or that immigrant populations could occupy ecological niches that developed following deglaciation and weathering of the rocks and minerals. At later times along the chronosequence, the turnover rate decreased to 0.006, indicating that the soil community composition had stabilized somewhat and the organisms present effectively prevented colonization by invading organisms. It is difficult to compare the turnover rates observed to those found in other studies of terrestrial habitats, because such measures are influenced by the sampling design, the area sampled and the metrics employed (Koleff et al. 2003). Some measures of turnover rate weigh the number of species present in both communities more, while others weigh more the number of species exclusive to one or the other location (Koleff et al. 2003). In addition, studies such as that of Noguez et al. (2005) employed a nested quadrat design, while we sampled along transects. Nonetheless, our results were generally consistent with those of Noguez et al. (2005) in which the turnover rate based on the species area relationship (Whittaker 1972; Arita & Rodriguez 2002) was determined to be 0.42 and 0.47 at two different locations in a biosphere reserve at the western coast in Mexico. These values were higher than those observed in animals, invertebrates and protists (Noguez et al. 2005).
The extraordinary diversity found in this glacial forefield is consistent with the results from other studies on bacterial diversity in polar environments and adds to the growing list of habitats, in which bacterial diversity has been shown to be unfathomable (Borneman et al. 1996; Torsvik et al. 1998; Zhou et al. 2004; Neufeld & Mohn 2005; Sogin et al. 2006; Elshahed et al. 2008). Additionally, similar to other studies (Roesch et al. 2007; ), significant fractions of the phylotypes in these communities were potentially novel. About 15% of all sequences could not be classified into known bacterial phyla and about 70% of all the sequences could not be classified into a known genus. Thus, with each additional study, we gain a better understanding of the full extent of prokaryotic diversity and its distribution on Earth. The study by Roesch et al. (2007) is the most comparable to ours, because the authors estimated bacterial richness in soils by analysis of 16S rRNA gene sequences obtained by pyrosequencing. Interestingly, Roesch et al. (2007) found that for a threshold of 97% sequence similarity phylotype richness in a Canadian forest soil (13 000 per g) was higher than that of three agricultural soils from temperate (Florida and Illinois, USA) and tropical (Brazil) climates, which were similar to one another (5–6 thousand per g). This was true for all thresholds of sequence similarity examined. Moreover, the species accumulation (rarefaction) curves for the Canadian soil markedly differed in shape from those determined for the other soils, and suggested there was greater phylotype evenness. Thus, the general finding of higher phylotype richness and evenness in the Canadian soil is reminiscent of that observed in the Arctic glacial forefield. Broad generalizations based on so few data are admittedly risky, but it is tempting to speculate that the latitudinal gradient in species diversity observed for plants and animals in which diversity is the greatest in the tropics and decreases towards the poles may not hold true for bacteria.
The difficulty faced in adequately sampling the high bacterial diversity in soils and differences in the computational methods used to estimate richness, evenness and species turnover rates hamper the assessment of prokaryotic species diversity and make it difficult to compare the results of different studies. For example, nonparametric methods used to estimate phylotype richness are known to be more conservative (Bunge 2008) and consistently result in lower phylotype richness estimates than parametric methods. This was the case in this study. Nonetheless, the estimates of phylotype richness both indicated that richness was high and followed the same trend over time, although the absolute values sometimes differed as much as two-fold (Fig. 1). It is not possible to know which estimates were closer to the true values, because the true underlying abundance distributions in the microbial communities were unknown. Similarly, different evenness indices may suggest different evenness of the members of the community, and in this study, we chose Evar because this index weighs the proportion of predominant members and rare members equally (Smith & Wilson 1996). Such discrepancies in estimates of diversity metrics are not uncommon because the means of estimating these metrics have limitations and assumptions that may not always be appropriate (for more detailed reviews see Beisel et al. 2003; Bunge & Barger 2008; Koleff et al. 2003; Smith & Wilson 1996). We have no easy solution for this problem and instead used multiple measures to obtain a range of possible values for each diversity metric.
As with other efforts to characterize bacterial diversity based solely on cataloguing 16S rRNA gene sequences, it is not possible to reliably infer the functions of the bacterial populations present. However, the general characteristics of some of the populations present in the glacial foreland soils are known and worth noting. Two phyla, Acidobacteria and Actinobacteria, were found in larger numbers along the chronosequences. Acidobacteria are commonly found in soils (Neufeld & Mohn 2005; Roesch et al. 2007) and are known to grow under oligotrophic conditions and at low pH (Dion 2008). Thus, the abundance of genera belonging to Acidobacteria in samples from the entire chronosequence and their increase in number over time may reflect the low levels of nutrients present in High Arctic soils and the decrease in soil pH that occurs along the chronosequence (Hodkinson et al. 2003). Several genera belonging to the phylum Actinobacteria were also found; in particular Nocardia found in our soil samples are known to form mycelia. The formation of mycelia is suggested to be advantageous in oligotrophic and desert environments, because the hyphae allow the bacteria to reach water and nutrients in pores further away (Dion 2008). More striking was the presence of N-fixing genera besides Cyanobacteria such as Rhizobium, Mesorhizobium and Bradorhizobium. Most studies in terrestrial polar environments on nitrogen fixation focused on Cyanobacteria (Whitton & Potts 2000); however, our data suggest that additional genera should be considered to contribute to N-fixation in High Arctic soils. In addition, a larger proportion of chemoorganotrophic phylotypes (32%) than photoautotrophs (10%) was observed at the youngest sites since glacier retreat, which contradicts the community assembly rule that photoautotrophic organisms are the first in early succession. These results are consistent with Skidmore et al. (2000), who detected viable chemoorganotrophic bacteria under the John Evans glacier, Canada and the suggestion of Hodkinson et al. (2002) that heterotrophic organisms may precede autotrophs in community assembly.
The mechanisms that work to maintain high bacterial diversity in soils of glacial forelands in the Arctic are unknown. We postulate that three mechanisms working individually or in some combination seem plausible. First, all soils are highly structured environments composed of soil aggregates that have a plethora of different microenvironments. These microenvironments differ from one another in many ways including pore water chemistry, the composition of particulates, the kind and amounts of organic matter, gas exchange, water holding capacity, the efficiency of binding bacterial cells and other features (Ehrlich & Newman 2008). Moreover, differences in micro-topography are associated with variation in temperature and wind. The matrix of pertinent environmental parameters creates a high order of ecological niches that can select for genotypes that are particularly well suited for local conditions. Second, bacterial populations in glacier foreland soils experience disturbances of intermediate strength on a range of temporal scales. For example, on a diurnal scale soil, temperatures can fluctuate from −10 to 25 °C within a few hours (Torsvik & Ovreas 2008), while over longer periods, there are changes in precipitation and in the glacial drainage system. The intermediate-disturbance hypothesis (Morin 1999) posits that such disturbances prevent the domination of a few organisms and the resulting community is even and species rich. Finally, competitive exclusion may be ineffective because of low bacterial growth rates that restrict the ability of one or more populations to exploit existing conditions and increase in number. Bacterial growth rates in High Arctic environments may be restricted because of lack of water, solar radiation, fluctuations in temperature and continuous freeze–thawing (Torsvik & Ovreas 2008). Moreover, temperatures are generally low and there is a paucity of nutrients that would support the growth of heterotrophs. The exceptionally slow growth of organisms in polar environments is illustrated by the fact that cyanobacterial cells in Antarctica were estimated to divide once every other year (Fritsen & Priscu 1996). Given the conditions that prevail, potential competitive advantages are counteracted by other factors that limit growth. Thus, individual populations are left to coexist and are unable to competitively exclude other phylotypes in the bacterial community.
Acknowledgements
The field work was funded by the Amundsen Center at the University of Tromsø and the Norwegian Polar Institute. The DNA Sequence Analysis Core Facility at the University of Idaho is supported by an NIH Center of Biomedical Research Excellence grant (P20 RR016448) from the National Center for Research Resources to LJF. Phylotype richness estimation was funded by NSF grant #0816638 ‘Estimating Biological Diversity and its Patterns’. This research was conducted using the resources of the Cornell University Center for Advanced Computing, which receives funding from Cornell University, New York State, the National Science Foundation, and other leading public agencies, foundations and corporations. This research was (partially) supported by Grant No. AB 2046 from the Student Grant Program at the University of Idaho, Moscow, ID. We thank Dr Rolf A. Olsen and the University Centre in Svalbard (UNIS) for facilitating our field research and Dr Ian Hodkinson and Dr Steve Coulson for information on their sampling locations and valuable discussions. We also wish to acknowledge Dr Eva Top for advice on methods for DNA isolation, Dr Chris Williams for valuable discussions on the experimental design, Dr Stefano Ventura and Silvia Turichia for their assistance in the field. Finally, we appreciate helpful comments on the manuscript provided by Jacob Pierson.
Conflicts of interest
The authors have no conflict of interest to declare and note that the sponsors of the issue had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.
