Research Article

Responses of tundra soil microbial communities to half a decade of experimental warming at two critical depths

Eric R. Johnston, Janet K. Hatt, Zhili He, Liyou Wu, Xue Guo, Yiqi Luo, Edward A. G. Schuur, James M. Tiedje, Jizhong Zhou, and Konstantinos T. Konstantinidis
PNAS July 23, 2019 116 (30) 15096-15105; first published July 8, 2019; https://doi.org/10.1073/pnas.1901307116

  1. Edited by Edward F. DeLong, University of Hawaii at Mānoa, Honolulu, HI, and approved June 10, 2019 (received for review January 23, 2019)

Significance

Ongoing permafrost thaw is expected to stimulate microbial release of greenhouse gases, threatening to further exacerbate climate change (cause positive feedback). In this study, a unique field warming experiment was conducted in Interior Alaska to promote surface permafrost degradation while maintaining uniform hydraulic conditions. After 5 winters of experimental warming by ∼1 °C, microbial community shifts were observed at the receded permafrost/active layer boundary, which reflected more reduced conditions, including increased methanogenesis. In contrast, increased carbohydrate utilization (respiration) was observed at the surface layer. These shifts were relatable to observed increases in CO2 and CH4 release from this study site and the surrounding ecosystem. Collectively, our results demonstrate that microbial responses to warming are rapid and identify potential biomarkers that could be important in modeling.

Abstract

Northern-latitude tundra soils harbor substantial carbon (C) stocks that are highly susceptible to microbial degradation with rising global temperatures. Understanding the magnitude and direction (e.g., C release or sequestration) of the microbial responses to warming is necessary to accurately model climate change. In this study, Alaskan tundra soils were subjected to experimental in situ warming by ∼1.1 °C above ambient temperature, and the microbial communities were evaluated using metagenomics after 4.5 years, at 2 depths: 15 to 25 cm (active layer at outset of the experiment) and 45 to 55 cm (transition zone at the permafrost/active layer boundary at the outset of the experiment). In contrast to small or insignificant shifts after 1.5 years of warming, 4.5 years of warming resulted in significant changes to the abundances of functional traits and the corresponding taxa relative to control plots (no warming), and microbial shifts differed qualitatively between the two soil depths. At 15 to 25 cm, increased abundances of carbohydrate utilization genes were observed that correlated with (increased) measured ecosystem carbon respiration. At the 45- to 55-cm layer, increased methanogenesis potential was observed, which corresponded with a 3-fold increase in abundance of a single archaeal clade of the Methanosarcinales order, increased annual thaw duration (45.3 vs. 79.3 days), and increased CH4 emissions. Collectively, these data demonstrate that the microbial responses to warming in tundra soil are rapid and markedly different between the 2 critical soil layers evaluated, and identify potential biomarkers for the corresponding microbial processes that could be important in modeling.

Representing only ∼16% of Earth’s terrestrial surface, northern-latitude permafrost soils and their overlaying active layers harbor an estimated 1,672 to 1,832 Pg of carbon (C), which accounts for ∼50% of the global soil organic C (SOC) reservoir (1, 2). This large C stock has accumulated and been preserved for thousands of years, primarily due to low temperatures and frozen conditions, which constrain microbial SOC mineralization (3, 4). Elevated atmospheric greenhouse gas concentrations are increasing global temperatures, and northern-latitude areas are experiencing a rate of warming that is more than twice the global average (5). As a result, regionally widespread and ongoing permafrost thaw is being observed (6⇓⇓–9). It has been estimated that permafrost could recede further by 30 to 70% by the end of the 21st century (10, 11). This alleviation of prior abiotic constraints is expected to stimulate microbial processes resulting in the release of greenhouse gases, primarily CO2 and CH4, which could further exacerbate climate warming (i.e., cause positive feedback) (12⇓–14). However, efforts to delineate microbial processes to improve future climate change predictions are hampered by the enormous complexity and heterogeneity of soil microbial communities and knowledge gaps regarding the microbial metabolic functions and their controls that operate in situ.

Furthermore, performing field warming experiments at northern latitudes is particularly challenging due to the remoteness of these sites and weather conditions. Previous investigations of laboratory-incubated permafrost soils under elevated temperatures reported shifts in community structure and functioning toward increased carbon respiration, even over a short period of a few weeks (15, 16). The relevance of these laboratory findings for in situ processes, however, remains speculative because the laboratory incubations cannot simulate closely the complexity of the natural environment, especially at the deeper soil layers that harbor older carbon stocks. Increasing temperatures and active layer thickness (i.e., an extension of maximum annual thaw depth) also affect other components of the soil, such as redox conditions, water availability, nutrient cycling, and aboveground vegetation (17⇓⇓–20). These physicochemical alterations, in turn, can result in additional, indirect effects of warming to soil microbes and their activities. However, such effects remain poorly understood and quantified to date.

More specifically, gradual shifts in vegetative communities can alter soil physical conditions, establish new plant–microbe symbioses, promote oxygen and or methane transit, and change the amount and type of plant material returned to soil (21⇓⇓–24). A combination of soil warming and redox conditions favorable for energy-yielding activities can result in a more rapid release of CO2 from tundra soils (25, 26). Under lower potential redox conditions, microbes might degrade SOC by gradually using electron acceptors with declining redox potentials, leading to conditions more suitable for methanogenesis. This can result in lower C loss compared with aerobic conditions and yet have greater feedback to climate due to larger CH4:CO2 ratio of emissions (26). Hence, the microbial and environmental features that serve as potential candidates for ecosystem-level evaluation of soil C release are vast and interconnected. Elucidating the key microbial processes and their environmental controls operating in situ requires an integrated approach that takes into account and measures the above parameters. Importantly, which of the microbial processes mentioned above—for example, aerobic/respiring vs. reducing/methanogenic pathways—prevail in northern latitudes and how fast these processes could be accelerated (or not) by warming remain essentially unknown. An improved understanding of the microbial clades and activities involved could help predict greenhouse gas fluxes and lead to an improved assessment of ecosystem-scale responses to climate warming (27).

To advance these issues, the Carbon in Permafrost Experimental Heating Research (CiPEHR) site was established in a discontinuous permafrost region in Interior Alaska in 2008 and represents the first in situ experiment to promote surface permafrost degradation without delaying spring soil melt (28, 29). In warmed plots, soil is insulated from frigid conditions during the winter months with snow that accumulated on the windward side of plastic mesh snow fences built by our team (SI Appendix, Fig. S1). This results in soil temperatures that are elevated year-round and also in the warming of deep soil, a manipulation not typically achieved with previous tundra warming experiments. Previous reports by our team demonstrated rapid physiological responses of belowground microorganisms after just 2 winters of experimentation at this site. However, community-wide shifts in microbial structure were less discernible at this early stage, perhaps due to a lag between physiological responses and corresponding shifts at the DNA level (30, 31). In this study, shotgun-metagenomic sequencing was used to assess responses of soil community structure after 4.5 y of experimental warming and gain insights into the undergoing in situ processes at 2 critical depths: 15 to 25 cm (preexisting active layer) and 45 to 55 cm (approximately the active layer/permafrost boundary at the initiation of experiment, which has since receded downward). The specific hypotheses tested with these data were: (i) microbial communities at the receded permafrost boundary layer depth will be more sensitive to elevated temperatures; (ii) at both soil depths studied, community functions stimulated by increased temperatures will include traits for the release of CO2 and/or CH4; and (iii) warming-induced changes in the relative abundances of traits involved in CO2 or CH4 release will be relatable to increased gas fluxes measured at the CiPEHR site and the surrounding region.

Results

Environmental Indices.

Experimental warming by means of increased winter snow cover resulted in a year-round elevation in soil temperature. For example, warming treatment plots were 1.02 and 1.11 °C warmer than control plots at 15- to 25-cm and 45- to 55-cm soil depths, respectively, during the growing season before sampling (from May to September; sampling took place in May; linear mixed-effects model [LME] used in conjunction with Tukey’s honest significant difference [HSD] test, P < 0.05) (Tables 1–3). Experimental warming also increased annual thaw duration from 104.2 to 112.2 d at the 15- to 25-cm depth (7.7% increase; nonsignificant with the LME model, but P < 0.05 using paired t test) and from 45.5 to 79.3 d in the 45- to 55-cm depth (74.4% increase; LME, P < 0.05). The mean thaw depth at the time of sampling in May 2013 was 18.3 and 23.0 cm in control and experimentally warmed plots, respectively (LME, P < 0.05). Experimental warming also increased mean aboveground plant biomass by 25.2%, increased volumetric water content slightly, and increased bulk density at the 45- to 55-cm depth profile (LME, P < 0.05). For aboveground plants, warming decreased lichen biomass by 41.8% on average (paired t test, P < 0.05), and increased the average biomass of Eriophorum vaginatum and Rhododendron tomentosum by 81.2% and 23.3%, respectively (paired t test, P < 0.05).

View this table:
Table 1.

Temperature and annual thaw duration of the experimental groups of the study

View this table:
Table 2.

Summary of soil physiochemical measurements for the experimental groups of the study

View this table:
Table 3.

Summary of soil physiochemical measurements for the experimental groups of the study

Broad Microbial Community Indices.

On average, 8.1 Gbp (billion base pairs) of sequencing effort per sample was obtained for 12 soil communities representing the 15- to 25-cm soil depth (6 warming, 6 control) and 5.8 Gbp per sample for 12 communities representing the 45- to 55-cm depth (SI Appendix, Table S1). Using Nonpareil 3 (options: -T kmer -k 32 -X 100000), the average estimated coverage of the microbial communities assessed with sequencing was 0.52 for 15- to 25-cm and 0.62 for 45- to 55-cm soil communities (SI Appendix, Fig. S2) (32). These coverage values suggest that, beyond the sequencing depth achieved, there is 52% and 62% likelihood that the additional sequences obtained would be redundant with ones already observed. The level of coverage obtained here (52 to 62%) should be appropriate for comparisons (e.g., low false-negative rate), with unsampled diversity presumably representing comparatively rarer taxa (33, 34). Experimental warming resulted in a near-significant increase in Nonpareil-derived sequence diversity Nd at the 45- to 55-cm depth, a metric that correlates well with α-diversity assessed by traditional metrics and is advantageous for shotgun metagenomic data (paired t test, P < 0.1) (SI Appendix, Fig. S3). MicrobeCensus was used to evaluate whether or not shifts in average genome size of the community due to the experimental warming could account for observed shifts in α-diversity or functional gene content (35). Such effects as a result of warming were not observed, but our analysis did reveal smaller genomes, on average, for microbial populations comprising the 45- to 55-cm soil depth compared with those at 15 to 25 cm (average of 4.7 Mbp vs. 6.7 Mbp; LME, P < 0.001) (SI Appendix, Table S1).

Soil Community Taxonomic Composition.

The relative 16S rRNA gene-based sequence abundance (metagenome derived) of Archaea was 1.69 ± 0.09% (mean ± SEM) and 2.64 ± 0.21% in 15- to 25-cm and 45- to 55-cm soil depths, respectively (LME, P < 0.05) (Fig. 1B and SI Appendix, Table S2). The archaeal community at the 15- to 25-cm depth was mostly represented by phyla Thaumarchaeota, followed by Euryarchaota at 1.11 ± 0.16% and 0.42 ± 0.13%, respectively. For the 45- to 55-cm layer, archaeal 16S rRNA gene sequences were mostly represented by Euryarchaeota and Crenarchaeota at 2.31 ± 0.22% and 0.23 ± 0.05%, respectively. Warming had no significant effect on the relative abundance of Archaea at the 15- to 25-cm depth (LME, P > 0.1), but increased archaeal abundance by a mean of 39% at 45 to 55 cm (2.21 ± 0.24% in control vs. 3.07 ± 0.27% in warmed soils; LME, P < 0.05). This increase was primarily driven by an increase in the abundance of the archaeal phylum Euryarchaeota (1.90 ± 0.27% vs. 2.73 ± 0.25%), and more specifically, the order Methanosarcinales, which was ∼3-fold more abundant in experimentally warmed plots at this soil depth (0.54 ± 0.14% vs. 1.61 ± 0.38%; LME, P < 0.05). Notably, a very strong correlation was observed between the relative abundance of Methanosarcinales and annual thaw duration (squared root regression R2 = 0.895, P < 0.001) (Fig. 1C). Most of the 18 bacterial phyla with mean relative abundances ≥0.1% across all 24 soil samples differed significantly between 15- to 25-cm and 45- to 55-cm soil communities in relative abundance (excluding only Thermotogae, Lentisphaerae, Armatimonadetes, and Actinobacteria; LME, P < 0.05) (Fig. 1A and SI Appendix, Fig. S4). At 15 to 25 cm, warming only significantly increased the relative abundance of Verrucomicrobia (8.0 ± 0.7% vs. 10.4 ± 0.6% relative abundance in control and warmed samples, respectively; P < 0.05). In contrast, at 45 to 55 cm, experimental warming resulted in significant shifts (P < 0.05) to phyla Proteobacteria (23.4 ± 2.4% vs. 30.4% ± 3.4), Actinobacteria (31.0 ± 4.1% vs. 20.6 ± 3.0%), and Microgenomates (0.2 ± 0.05% vs. 0.8 ± 0.2%).

Fig. 1.
Fig. 1.

Taxonomic shifts as an effect of experimental warming. (A) Mean relative abundance of bacterial phyla for each depth × treatment combination. Underlying data are based on 16S rRNA gene-encoding fragments recovered from metagenomic datasets. Values represent the abundance of each bacterial phylum as a proportion of the total bacterial community. Only phyla with a mean relative abundance higher than 0.1% across all 24 datasets are displayed. (B) Mean relative abundances of archaeal phyla for each depth × treatment combination. Underlying data are based on 16S rRNA gene-encoding fragments recovered from metagenomic datasets. Values represent the abundance (as a percentage) of each archaeal phylum relative to the total abundance of all recovered bacterial and archaeal 16S rRNA gene fragments (i.e., the total prokaryotic community). Error bars represent the mean ± the SEM (n = 6) for cumulative (total) relative Archaea abundance. (C) Correlation between the relative abundance of Methanosarcinales and annual thaw duration (in days of the year) for 45- to 55-cm soils. Linear regression was fitted for n = 12 points. Significance of the correlation coefficient was determined using a two-tailed Student’s t-distribution (r = 0.946, df = 11).

Consistent with the abovementioned results, at 45 to 55 cm, experimental warming resulted in significant community-wide shifts in taxonomic β-diversity based on abundance weighted Jaccard distances of phylum-level community composition (analysis of similarities [ANOSIM] and multiresponse permutation procedures [MRPP], P < 0.05) (Fig. 2 and Table 4). There was also a significant relationship between β-diversity (derived from either UniFrac or Jaccard distances) and annual thaw duration or depth of thaw at time of sampling (Adonis, P < 0.05). There were no significant relationships between experimental treatment or thaw indices and 16S rRNA gene-based or phyla-based evaluations of community structure for samples representing the 15- to 25-cm depth, similar to observations of this depth made after just 2 winters of experimental warming (30).

Fig. 2.
Fig. 2.

Functional and phylogenetic shifts as an effect of experimental warming. (A) Principle coordinates analysis (PCoA) plot of consolidated KO term annotations. Underlying data are based on abundance-weighted jaccard distance matrix derived from a KO term counts matrix. (B) PCoA plot of community phylogenetic composition. Underlying data are a weighted unifrac distance matrix of 16S rRNA gene-encoding fragments recovered with Parallel-META and processed in the QIIME software package, as described in Materials and Methods. (C) PCA plot of Swiss-Prot gene annotations. Underlying data are based on a gene count matrix consolidated from Swiss-Prot database references, which underwent variance-stabilizing transformation using the DESeq2 package. The 24 samples shown in the 2D plane of each plot are spanned by their first 2 principal components.

View this table:
Table 4.

The effects of warming and thaw depth or duration on soil microbial community functional and phylogenetic structure

Shifts in Energy-Yielding Metabolisms with Experimental Warming.

In addition to the observed shifts in microbial community phylogenetic structure for 45- to 55-cm metagenomes, experimental warming shifted community functional gene content (summarized as Carbohydrate-Active enZYmes [CAZy] families or Kyoto Encyclopedia of Genes and Genomes, or KEGG, Orthology [KO] terms) (ANOSIM, P < 0.05 or P < 0.1) (Fig. 2 and Table 4). Also for this depth, β-diversity based on KO terms was significantly associated with annual thaw duration and measured thaw depth at time of sampling (Adonis, P < 0.05). For 15- to 25-cm soils, warming only significantly altered community β-diversity reflecting the composition of genes involved in carbohydrate metabolism and binding (i.e., CAZy protein families; MRPP, P < 0.05; ANOSIM, P < 0.1), and no significant relationships between functional structure and thaw indices for this depth were identified.

Furthermore, for 15- to 25-cm communities, experimental warming increased the relative proportion of functional genes matching to CAZy reference sequences by 7.9% (LME, P < 0.05) (Fig. 3A). This included a 12.0% increase in glycoside hydrolase genes, a 2.6% increase in glucosyltransferase genes, an 8.1% increase in carbohydrate binding genes, an 8.9% increase in carbohydrate esterase gene, and a 20.8% increase in polysaccharide lyase genes (LME, P < 0.05) (Fig. 3 B and C). These results were generally consistent with Gene Ontology (GO) process categories involving the catabolism of various organic matter substrates (Fig. 3D). Of the 25 categories based on these definitions that differed by ≥10% between control and warmed 15- to 25-cm soils, 22 of 25 were more abundant in experimentally warmed metagenomes, with 4 having an adjusted P < 0.05 with LME. Similar shifts were not found between warmed and control 15- to 25-cm soils from the 1.5-y collection (SI Appendix, Fig. S5). While warming also shifted community β-diversity of CAZy families at the 45- to 55-cm depth, there was only a significant 3.1% increase in the abundance of glucosyltransferase genes at this depth (LME, P < 0.05).

Fig. 3.
Fig. 3.

Shifts in carbohydrate metabolism genes as an effect of experimental warming. Relative abundances of (A) all CAZy annotations, and CAZy modules (B) glycoside hydrolases (GH), glycosyl transferases (GT), and (C) auxiliary activities (AA), carbohydrate binding (CBM), carbohydrate esterases (CE), and polysaccharide lyases (PL). Underlying values represent the mean relative abundance of each category for each sample group (i.e., treatment × depth). The abundances of CAZy functions were determined by normalizing the number of annotations matching to each broad definition by the number of annotations matching to the more comprehensive Swiss-Prot database. Letters distinguish sample groups that were significantly different (adjusted P < 0.05); that is, values with letters differing from letters assigned to other groups designate a statistically significant difference between groups. Significance was determined using a LME model where experimental fence was treated as a within-subjects factor, in conjunction with Tukey’s HSD test. (D) Shifts in the relative abundances of genes involved in the catabolism of various organic matter substrates with 4.5 y of experimental warming at 15 to 25 cm. The relative abundances of catabolic pathways were determined by matching Swiss-Prot references to GO Biological Process terms, and dividing the number of annotations for each term by the total number of Swiss-Prot annotations. Bars reflect the mean relative abundance for each sample group (unwarmed, warmed) at 15 to 25 cm (n = 6) and error bars represent the SEM. Only pathways that differed between warmed and unwarmed soils by ≥10%, on average, are shown.

Warming of the 15- to 25-cm soil layer did not result in any KO term that changed in abundance with an adjusted P < 0.05 (using DESeq2) (Datasets S1 and S2). However, of the 3,274 KO terms with abundances adequate for P value assignment in DESeq2, 162 KO terms differed significantly (adjusted P < 0.05) and another 153 differed near-significantly (adjusted P < 0.1) between control and warmed soil community samples at 45 to 55 cm. Also for the 45- to 55-cm depth, the abundances of 490 KO terms increased or decreased significantly across a gradient of annual thaw duration, and another 301 KO terms displayed a near-significant change with thaw duration. These shifts included a broad increase in functions involved in methanogenesis (or potentially, reverse methanogenesis) (36) between control and experimentally warmed communities (Fig. 4A). This evaluation was primarily based on KEGG modules M00567, M00357, M00356, and M00563, corresponding to methanogenesis from CO2, acetate, methanol, and methylamines (37) (Fig. 4A). Among these shifts, the largest and consistently significant responses were observed for KO terms specific to methanogenesis from acetate, which increased on average by 170%. There were also strong correlations between the abundance of archaeal order Methanosarcinales and KO terms for methanogenesis from acetate (R2 = 0.93; P < 0.01).

Fig. 4.
Fig. 4.

Shifts in microbial energy-generating pathways as an effect of experimental warming. (A) Heatmap showing the relative abundances of genes involved in methanogenesis (or potentially, reverse methanogenesis) for 45- to 55-cm depth soil metagenomes. KO terms under the “Core methanogenesis genes” category represent those used by all methanogens. Other categories and KO terms refer to pathways involved in the usage of different electron acceptors to sustain methanogenesis, which are variable between methanogenic taxa. To emphasize differences between datasets, the values for each gene were normalized by the mean relative abundance of all 12 of the 45- to 55-cm soil metagenomes (by rows). Subplots on the right of the heatmap represent the log2 fold-difference calculated for each gene, and are colored by whether differences between control and experimentally warmed soils were statitistical significant (Right) or if there was a significant association between gene relative abundance and annual thaw duration (Left) (dark gray = adjusted P < 0.05; light gray = nonadjusted P < 0.05; white = nonsignificant). Subplot on the top of the heatmap represents the annual thaw duration (i.e., days of the year soil was thawed) for each plot represented by the corresponding metagenomes. (B) Heatmaps showing the relative gene abundances for assimilatory and dissimilatory nitrate and sulfate reduction. Underlying data for A and B are communities respresenting 45- to 55-cm soil samples.

Increased temperatures or extended annual thaw of 45- to 55-cm soils also increased the relative abundances of genes involved in dissimilatory sulfate reduction and dissimilatory nitrate reduction to ammonium (adjusted P < 0.05 or 0.1 with thaw or treatment) (Fig. 4B), but resulted in no significant changes in denitrification genes norB and nosZ. The opposite trends were observed for KO terms specific to assimilatory nitrate reduction and assimilatory sulfate reduction (Fig. 4B), implying more reducing redox conditions upon warming and permafrost thawing. An increase in the abundances of genes involved in dissimilatory sulfate reduction was largely attributable to the Proteobacteria order Syntrophobacterales, a clade with many known sulfur-reducing taxa, which increased in abundance from 1.0 to 2.5% with experimental warming (P < 0.005). Consistent with this, a recovered metagenome-assembled genome (MAG) belonging to Syntrophobacterales was found to possess genes for respiratory sulfur reduction (dsrAB, aprAB).

DESeq2 analysis revealed statistically significant decreases in the relative abundances of cytochrome c oxidase genes with an average 2012 growing season temperature for 45- to 55-cm soils, even when controlling for thaw indices, consistent with more reducing conditions. Multiple linear regression analysis revealed that average summer 2012 temperature and thaw depth had significant independent (controlling for the other variable) as well as additive negative associations with the abundances of cytochrome c oxidase genes (multiple R2 > 0.8; P < 0.001) (Fig. 5A), and conversely, positive associations with genes involved in methane production from acetate (multiple R2 = 0.85; P < 0.0005) (Fig. 5B). The total relative abundance of all CAZy functions in 15- to 25-cm communities, as well as those belonging specifically to carbohydrate binding or glycoside hydrolase modules, had significant positive correlations with cumulative ecosystem respiration in the month following sampling (June 2013) (t-distribution, r ≥ 0.6, P < 0.05) (see also Fig. 5C).

Fig. 5.
Fig. 5.

Shifts in energy generating metabolisms with temperature and thaw depth in 45- to 55-cm soil communities. (A) Independent and additive effects of thaw depth and average summer temperature on the relative abundances of cytochrome c oxidase genes and (B) genes involved in methane generation from acetate. (C) Cumulative June 2013 ecosystem respiration for each experimental plot. Three-dimensional plots were made using the R package “Scatterplot3d” and illustrate relationships between average summer soil temperatures, thaw depth at time of sampling, and the abundance of energy generating metabolisms summarized as KO terms (the title of each plot reflects the KO identifier, as well as the common gene identifier). R function lm() was used to fit multiple variables to determine assess their independent and combinatorial effects. Individual data points are colored by average summer 2012 temperature, which proceed from blue to red to reflect low vs. comparatively higher soil temperatures.

Recovery of MAGs.

Assembly and population genome binning led to the recovery of 173 MAGs with a quality score ≥60 (calculated as completeness – 5 × contamination based on CheckM) (38, 39) (SI Appendix, Tables S3 and S4). These medium- to high-quality MAGs collectively recruited 17.0% and 24.9% of the short reads, on average, for the 15- to 25-cm and the 45- to 55-cm depth, respectively (using default megablast alignment; ≥97% nucleotide ID and ≥100-bp alignment). Dereplication of highly similar MAGs based on FastANI (40) values of ≥95% resulted in the consolidation of 74 nonredundant MAGs. This included 4 nonredundant archaeal MAGs, all derived from 45- to 55-cm soil metagenomes (SI Appendix, Fig. S6). Two archaeal MAGs were identified as belonging to order Methanosarcinales and the other two were assigned to Methanocellales based on taxonomic identification with the Microbial Genomes Atlas (MiGA) (SI Appendix, Fig. S7). Each archaeal MAG possessed genes involved in methanogenesis using CO2, acetate, and in one case, methanol, as electron acceptors. One Methanosarcinales MAG matched closely to ANME-2d taxa Candidatus Methanoperedens sp. BLZ1 and BLZ2 (87.5% average nucleotide identity [ANI]) (41, 42) and was also assigned to genus Methanoperedens with GTDB-Tk (SI Appendix, Table S4). Similar to the previously recovered Candidatus Methanoperedens taxa, the CiPEHR MAG possessed genes for reverse methanogenesis, nitrate reduction to nitrite (narGH), a cytochrome c nitrite reductase (nrfA), as well as nitrogenase genes nifHDK.

Among bacterial MAGs that shifted in abundance with warming, there was a variable yet large 79% increase in the most dominant bacterial population (making up ∼1–2% of the total community) in the 15- to 25-cm active layer with warming (paired t test, P < 0.1), which represented a member of the Acidobacteriaceae family (53.1% amino acid identity [AAI] with Candidatus Koribacter versatilis Ellin345 as determined with MiGA; assigned to family Koribacteraceae with GTDB-Tk) (SI Appendix, Fig. S8 and Table S4). These results were consistent with smaller, nonsignificant shifts observed after 1.5 y of experimentation for the same taxon. This dominant MAG was previously found to be widespread throughout Alaskan tundra and appears to encode diverse metabolisms for labile and recalcitrant organic matter degradation (31).

Discussion

Temperature increases as a result of experimental warming lasted year-round and were moderately uniform between the depth profiles evaluated here (15 to 25 and 45 to 55 cm) (Tables 1 and 2). Climate change and increased temperatures will undoubtedly change numerous important aspects of tundra ecology, such as thaw depth and plant community factors, which can lead to indirect effects of warming on soil microbial community functioning (18, 25, 43). The magnitude and direction of the microbial responses to these factors could vary considerably between depth profiles. Thus, warming of deep soil extending down to the permafrost boundary layer was a desirable outcome of our experimental warming manipulation. Experimental warming also increased annual thaw duration by 7.7% (+8 d) (relative to control plots) at the 15- to 25-cm soil layer and, in strong contrast, by 74.4% (+33.8 d) at the 45- to 55-cm layer for the year preceding sample collection.

Consistently, greater overall changes in community structure were observed for 45- to 55-cm soil communities. This included an ∼3-fold increase in Methanosarcinales, a methanogenic order and the most abundant archaeal clade observed in these soils. The abundance of this Methanosarcinales was highly relatable to annual thaw depth (R2 = ∼0.9, P < 0.001) (Fig. 1C) and accompanied a comprehensive increase in the relative abundances of methanogenesis genes (Fig. 4A). The increase in the relative abundance of these methanogens corresponded to increased CH4 emissions at the CiPEHR site due to our experimental warming and nearby regions due to ambient warming. For example, CH4 release was considerably greater from experimentally warmed plots (44) in the same year that soil cores were collected for the current study, and the warming effect from CH4 emissions is now larger than that from CO2 in this general area (the Eight Mile Lake region) (45).

Aerobic methane oxidation genes were fairly low in relative abundance, and those that were over the detection limit of our metagenomics effort did not shift in response to warming at either depth (SI Appendix, Fig. S9). It is possible that methane was consumed at a depth not evaluated in this study or by anaerobic methane oxidizing microbiota that remain poorly understood (36, 46). Related to the latter, one of the recovered archaeal MAGs was identified as a close relative to recently described ANME-2d taxa (87.5% ANI to Candidatus Methanoperedens sp. BLZ1) (41), and possessed genes for reverse methanogenesis, dissimilatory nitrate reduction to ammonia (DNRA), and nitrogen fixation (nitrogenase genes nifHDK). Previous reports have demonstrated nitrogen fixation in ANME-2 Archaea (47, 48), but the presence of functions for both DNRA and N2-fixation in an ANME-2d genome is unique. This combination of traits (producing ammonium as a waste product of DNRA while also capable of fixation of gaseous N2 to ammonium, which is energetically expensive) could signify a high degree of metabolic versatility or alternative functions associated with these genes. Related to this, previous research has provided some evidence of ANME taxa capable of both forward and reverse methanogenesis (49); hence, it remains unclear what role the ANME-2d population has in methane fluxes in these tundra soils.

For the last decade (2008–2017), the majority of tundra habitats have exhibited a 2 to 4.2 °C temperature increase above the 1950–1980 base period during cold seasons of the Northern Hemisphere (November to April) (50, 51). Despite an increase in background temperature that has continued to recede the permafrost boundary layer in warmed and unwarmed plots, a rapid and clear response of methanogenic taxa was nonetheless observed in experimentally warmed plots. Elevated temperatures and increased thaw were associated with independent as well as additive declines in the abundances of O2-specific respiratory metabolisms, such as cytochrome c oxidase genes (Fig. 5A). The opposite trend was found for genes involved in methane production from acetate (Fig. 5B), reflecting a shift toward lower reduction potential. It is possible that accelerated warming caused by our experimental set up caused a more thorough depletion of electron acceptors compared with gradual warming (i.e., recent ambient), due to a smaller period over, which the acceptors can be replenished, resulting in a greater CH4:CO2 ratio of emissions. A decline in redox conditions was also evidenced by greater relative abundances of DNRA and dissimilatory sulfate reduction genes (Fig. 4B), and associated taxa, such as Proteobacteria order Syntrophobacterales. Accompanying declines in assimilatory nitrogen and sulfur reduction genes observed under experimentally warmed conditions at 45 to 55 cm could be due to enhanced NH4+ and S2− production from dissimilatory reduction pathways, and thus a selection of traits for direct acquisition of these reduced compounds (Fig. 4B).

Warming also increased the biomass of E. vaginatum by 81%, a plant species that excretes organic acids such as acetate that can fuel CH4 production, transit, and release (23). Another plant species known to fuel methanogens, Carex rostrata (23), is common in this field location, albeit with lower biomass. It is possible that a greater response observed for genes specific to methanogenesis from acetate compared with other methanogenesis pathways (Fig. 4A) is attributable at least in part to these types of plant–microbe interactions. However, the associations between plant and methanogens resulting in greater CH4 emission from tundra habitats remains elusive (52, 53), and thus are topics worthy of consideration in future studies.

At 15 to 25 cm, warming increased the total proportion of community functional genes involved in carbohydrate metabolism (Fig. 3A), including glycoside hydrolases, glycosyltransferases, carbohydrate esterases, carbohydrate binding, and polysaccharide lyases (Fig. 3 B and C). Furthermore, the relative abundances of these carbohydrate metabolism categories were relatable to increased ecosystem respiration as a result of experimental warming in the month following sampling (June 2013) (Fig. 5C) (r ≥ 0.6, P < 0.05). An early assessment of this depth profile after 1.5 y of experimental warming found a similar response using GeoChip that was undetected or less obvious with metagenomic analysis (30). Thus, the response observed here is consistent with earlier observations, but also implies an ongoing change over time (SI Appendix, Fig. S5). Consolidation of CAZy reference sequences into family definitions also revealed a significant relationship between warming treatment and the relative composition (i.e., β-diversity) of these functions (Table 4). Moderately strong associations between the fraction of genes involved in carbohydrate metabolism and the average 15- to 25-cm soil temperature during the same month for which soils were collected (May) (R2 = 0.521; P < 0.01) and the average temperature of the preceding winter season (R2 = 0.478; P < 0.05) were also observed. These direct associations between the community fraction of genes involved in carbohydrate metabolism and soil temperature were much stronger than correlations from other broad environmental measure displayed in Tables 1–3, including plant biomass (R2 = 0.179 or less; P > 0.1), which increased by 25.2% in experimentally warmed plots (Table 3) (paired t test P < 0.05). These results of stimulated functions involved in organic C turnover as a direct response to warming are somewhat surprising, as recent C releases in an adjacent area were found to result mostly from new SOC sources (17). This could suggest that even with enhanced vegetative production, increases in the relative abundances of genes involved in carbohydrate metabolism are nonetheless more attributable to soil temperature, at least in frigid locations where microbial catabolic functions are constrained by cold conditions.

Conclusions and Future Perspectives

Our evaluation revealed a much greater susceptibility of soil communities at the recently receded permafrost/active layer boundary to further thaw and temperature increase; these responses contrasted with warming-induced changes to shallower, preexisting active layer communities. Community metabolic shifts in the deeper layer strongly indicated a stimulation of methanogens, particularly those using acetate for CH4 production, which was consistent with independent measurements of methane emissions from the same site (44). These effects were not solely attributable to increased thaw depth or duration, but were in conjunction with increased temperatures. If accelerated warming results in more thorough depletion of electron acceptors than would occur under a more gradual thaw, it could result in greater CH4:CO2 of emissions, and thus stronger positive feedback to climate change. Meanwhile, increasing CH4 emissions observed in this area now outweigh the warming potential of CO2 emissions and offsets C uptake by local plants (44, 45). Furthermore, northern-latitude areas have experienced a rate of warming over the past several decades that is much greater than the global average. If responses under rapid vs. gradual warming are dissimilar in potential CH4 release, this could imply that inferences made from natural thaw gradients (54) or geological records may not serve as adequate predictors of future CH4 release from tundra soils if warming continues at a rapid pace. Shifts observed at the 15- to 25-cm layer were distinctly attributable to functions involved in carbohydrate metabolisms, were most attributable to the direct effect of elevated temperatures, and were also relatable to enhanced ecosystem respiration.

While an increase in the abundance of methanogenic taxa and genes certainly implies that these CH4-generating taxa were active at some point, numerous environmental conditions and biological factors (i.e., plants, methanotrophs) ultimately determine methane emission vs. consumption. One archaeal MAG closely matching the recently described ANME-2d taxa possessed genes for reverse methanogenesis, DNRA, and nitrogen fixation. This combination of traits highlights uncertainty regarding its role in N and C flux in tundra systems. Hence, there remains a need to understand how susceptible these soils are to CH4 loss and the microbial mechanisms involved. Assessments of in situ microbial activity, such as through the use of proteomics and transcriptomics of soils (16, 55, 56), can offer such insights regarding whole-community or individual-taxon functioning. The emerging breadth of approaches enabling a more comprehensive investigation of soil microbiota should further serve to unravel the complex responses and interactions between ecological attributes in warming tundra ecosystems. An improved understanding of warming tundra ecosystems will greatly improve our understanding of their contributions to ongoing climate change. The gene and genome sequences reported here should help facilitate proteomics and primer design for PCR assays that can be used to precisely monitor the dynamics of responsive populations in tundra and elsewhere, and further corroborate the results reported herein.

Materials and Methods

Study Site Description and Sample Collection.

The CiPEHR site was established at a moist acidic tundra area in September 2008 in Interior Alaska near Denali National Park in the Eight Mile Lake region (63°52′59″N, 149°13′32″W) in a discontinuous permafrost region where permafrost thaw has been observed in the past several decades. Experimental design and site description have been described in detail previously (28); see also SI Appendix, Fig. S1 for a diagram and photos of the experimental fences. In brief, 3 experimental blocks were located ∼100 m away from each other. In each block, 2 snow fences were constructed about 5-m apart in the winter. Snow fences were made out of plastic mesh netting with ∼2-inch2 holes. The winter warming treatment plots were located 5-m back from the leeward side of the snow fences, while plots at the windward side of snow fences served as paired controls. Elevated temperatures resulted from thicker snow cover on the soil surface and lower wind strength during the winter months. Accumulated snow was removed in the spring before snow melt to maintain uniform hydraulic conditions between winter warming and control treatments.

Air temperature from 2004 to 2013 ranged from a monthly average of −18.0 ± 1.8 °C in January to 13.4 ± 0.5 °C in July, with an average annual temperature of −2.7 ± 0.4 °C. The mean annual precipitation was 378 mm and the mean growing season precipitation from 2004 to 2013 was 216 ± 24 mm (44). Only C3 plant species were observed in this area; dominant species include E. vaginatum, Vaccinium uliginosum, some other vascular species, nonvascular feather moss, and lichen. In the experimental plots, the upper 45 to 65 cm of soil was rich in organic C materials and below was mineral soil with a mixture of glacial till and windblown loess. The active layer depth was about 50 cm at initiation of the experiment but has since expanded to lower depths (44, 45). Twelve soil cores, 6 each from warming treatment and control plots, were collected using electric drills in the beginning of the 2010 and 2013 growing season (May), 2 (1.5 y) and 5 winters (4.5 y) after the start of the winter warming manipulation. Cumulative ecosystem respiration for the month following sampling was based on half-hourly chamber measurements used to evaluate net ecosystem exchange, where nighttime net ecosystem exchange was taken as ecosystem respiration (57). Additional information on site monitoring and characterization of environmental indices can be found in SI Appendix.

Soil DNA Sequencing and Read Processing.

Soil DNA extraction, metagenome library preparation, and sequencing of 4.5-y soil samples were performed similarly to methods used for soils collected after 1.5 y (31); see SI Appendix for details. The 4.5-y metagenome datasets were deposited in the European Nucleotide Archive under study no. PRJEB31848 (58). Accession IDs corresponding to each sample metagenome are provided in SI Appendix, Table S1. Metagenomes representing the 1.5-y sample collection are available under ENA project ID PRJEB10725 (31). Metagenomic paired-end reads were merged using PEAR (59) (options: -p 0.001). All merged and nonmerged reads were then quality-trimmed with the SolexaQA package (60) (options: -h 17; ≥98% accuracy per nucleotide position). Trimmed sequences used downstream for functional annotation or taxonomic assignment were truncated to 150 bp to avoid read-length biases.

Soil Community 16S rRNA Gene Analysis.

For the assessment of taxonomic composition, 16S and 23S rRNA gene fragments were first recovered from metagenomes using SortMeRNA (61). The relative proportion of sequences matching to bacterial vs. archaeal rRNA gene sequences from SortMeRNA was used to estimate the relative abundances of either domain. The relative abundances of bacterial phyla were determined with Parallel-META 2.0 (62), using SortMeRNA output sequences as input. Archaeal phyla abundances were determined by aligning SortMeRNA output sequences against the SILVA SSU database (v132) (63) using BLASTN alignment (BLAST+ v2.2.28; options: -word_size 18) (64). Matches to archaeal sequences with ≥80% nucleotide identity and ≥100-bp alignment length were retained. To assess community phylogenetic structure, SortMeRNA output sequences were used with the QIIME 1 software package (65) for closed-reference OTU picking with UCLUST (options: -m uclust_ref -C -z) (66) against the SILVA SSU database dereplicated at 99% nucleotide identity. Aligned SILVA reference sequences represented in the OTU table were then used to construct a phylogenetic tree using FastTree (67). The resulting phylogenetic tree and OTU count table was then used to generate a weighed UniFrac distance matrix summarizing pairwise phylogenetic distances between communities.

Soil Community Functional Gene Assessment.

Short nucleotide sequences were searched against the Swiss-Prot (downloaded on November 27, –2016) (68) and CAZy (downloaded on July 15, 2016) (69) reference databases using DIAMOND BLASTX alignment (options: -k 1 -e 1E-5–sensitive) (70). Matches to Swiss-Prot or CAZy reference sequences with a bit score ≥55 were retained for further analysis. Independent count matrices were made to summarize the results of alignment against either database. A reference data file associating KO functions (KO terms) (37) to Swiss-Prot database entries obtained from https://www.uniprot.org/downloads was used to convert Swiss-Prot annotation counts into a more consolidated count matrix of KO terms (available as Dataset S1). Swiss-Prot annotations were also consolidated into GO categories as described previously (31). Annotations to CAZy reference sequences were consolidated into a count matrix of structurally related catalytic and carbohydrate-binding families. They were also consolidated further into the 6 broad functional modules covered in the CAZy database, including glycoside hydrolases, glycosyl transferases, auxiliary activities, carbohydrate binding, carbohydrate esterases, and polysaccharide lyases. For each sample metagenome, the abundances of summarized CAZy definitions, as proportions of all functional genes, were determined by normalizing the number of annotations matching to each broad definition by the number of annotations matching to the more functionally comprehensive Swiss-Prot database as well as the total number of sequences originally used as query for functional annotation. For the summarized count matrices of KO terms or CAZy enzyme families, compositional dissimilarity between samples was calculated with Bray–Curtis and abundance-weighted Jaccard distance metrics.

DNA Assembly, Genome Binning, and MAG Quality Assessment.

Metagenomic sequences were assembled with IDBA-UD (71) (options:–mink 43–maxk 123–step 4–min_contig 300). Resulting contigs ≥ 2.5 kbp were used with MetaBAT2 (options:–minCVSum 10) (72) and MaxBin2 (options: -min_contig_length 2000) (73) to recover microbial population genomes. Before binning, Bowtie 2 was used to align short-read sequences to embled contigs (options:–very-fast) (74) and SAMtools was used to sort and convert SAM files to BAM format (75). Sorted BAM files were then used to calculate the coverage (mean representation) of each contig in each sample metagenome. The quality of each resulting MAG was evaluated with CheckM v1.0.3 taxonomy workflow for Bacteria and Archaea, separately (39). The result from either evaluation (i.e., taxonomy workflow for Archaea or Bacteria) with the highest estimated completeness was retained for each MAG. MAGs with a quality score ≥60 were retained (from ref. 38; calculated as the estimated completeness – 5 × contamination). FastANI (40) was used to estimate the ANI between MaxBin 2 and MetaBAT 2 generated MAGs. For redundant MAGs obtained independently from either binning method that matched at ≥95% ANI and were derived from the same assembled metagenome, the MAG version with the highest quality score or greatest overall size was retained while the other was discarded. Protein-encoding genes from MAG contig sequences were predicted with Prodigal (76) and the resulting nucleotide sequences were searched against the Swiss-Prot database reference sequences using DIAMOND BLASTX alignment as described above. Select MAGs were also processed using the MiGA webserver for taxonomic assignment (77). MAGs uploaded to the MiGA online webserver are freely available under MAGs project “AK Tundra” for browsing of their amino acid identity relatedness, taxonomic classification against named taxa, and quality of genome sequence (http://enve-omics.ce.gatech.edu:3000/projects/22). MAG taxonomic assignments with MiGA and European Nucleotide Archive accession IDs are provided in SI Appendix, Table S5. GTDB-Tk was also used for taxonomic assignment of MAGs (SI Appendix, Table S4) and for generating phylogenetic trees for bacterial (Dataset S3) and archaeal (Dataset S4) MAGs (78).

Statistical Analyses.

With the R package nlme (79), a LME model treating experimental fence as a within-subjects factor was used to evaluate differences between experimental groups (15- to 25-cm unwarmed, 15- to 25-cm warmed, 45- to 55-cm unwarmed, 45- to 55-cm warmed) based on their environmental indices, the relative abundances of Bacteria, Archaea, and their phyla, and the relative proportions of summarized CAZy modules or broad GO biological process terms. Tukey’s HSD test was used as a post hoc analysis for assessing statistical significance of differences in the means of 4 treatments (15- to 25-cm control, 15- to 25-cm warming, 45- to 55-cm control, and 45- to 55-cm warming). ANOSIM and MRPP were used to evaluate the significance of compositional differences between warmed and control community functional and phylogenetic structure at either depth. Adonis was used to identify statistically significant associations between sample distances and pertinent environmental indices (e.g., thaw depth at sampling, annual thaw duration). Significance of correlations between metadata and certain community measures (metagenome annotations) was assessed using a 2-tailed t-distribution. Tables with raw counts of KO terms were processed with the DESeq2 software package (80) to identify significant, differentially abundant functions between treatments within a depth (where fence was treated as a within-subjects factor using design = ∼ Fence + Treatment) or using annual thaw duration as a continuous variable. Benjamini–Hochberg P value correction (81) was used to account for false-discovery rate arising from multiple comparisons. DESeq2 results are provided as Dataset S2. The R package Superheat (82) was used to illustrate the relative abundances of annotations in conjunction with metadata and statistical results. The lm() function in R was used to perform multiple linear regression to evaluate potential interactive effects between thaw depth and average summer 2012 temperature on the relative abundances of genes involved in energy-yielding metabolisms. Paired t test was used to compare certain ecological factors that did not involve a depth component.

Acknowledgments

This research was supported by US Department of Energy Award DE-SC0004601 (to J.Z., Z.H., L.W., Y.L., E.A.G.S., J.M.T., and K.T.K.) and US National Science Foundation Awards 1356288 and 1759831 (to K.T.K.). We thank Joel Kostka for his feedback on this manuscript as well as the Partnership for an Advanced Computing Environment at the Georgia Institute of Technology, which enabled the computational tasks associated with this study.

Footnotes

  • ↵1To whom correspondence may be addressed. Email: kostas{at}ce.gatech.edu.

  • Author contributions: E.R.J., Z.H., L.W., Y.L., E.A.G.S., J.M.T., J.Z., and K.T.K. designed research; E.R.J., J.K.H., and X.G. performed research; E.R.J. analyzed data; and E.R.J. and K.T.K. wrote the paper.

  • The authors declare no conflict of interest.

  • This article is a PNAS Direct Submission.

  • Data deposition: Metagenome and MAG sequences have been deposited in the European Nucleotide Archive (study ID PRJEB31848). Individual accession IDs for metagenomes and metagenome-assembled genomes can be found in SI Appendix, Tables S1 and S5, respectively.

  • This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1901307116/-/DCSupplemental.

Published under the PNAS license.

View Abstract
Back to top
Article Alerts
Email Article
Citation Tools
Request Permissions
Responses of tundra soil microbial communities to half a decade of experimental warming at two critical depths
Eric R. Johnston, Janet K. Hatt, Zhili He, Liyou Wu, Xue Guo, Yiqi Luo, Edward A. G. Schuur, James M. Tiedje, Jizhong Zhou, Konstantinos T. Konstantinidis
Proceedings of the National Academy of Sciences Jul 2019, 116 (30) 15096-15105; DOI: 10.1073/pnas.1901307116
Digg logo Reddit logo Twitter logo Facebook logo Google logo Mendeley logo
Proceedings of the National Academy of Sciences: 116 (30)
Table of Contents

Submit

Article Classifications

Jump to section

You May Also be Interested in

Large piece of gold
News Feature: Tracing gold's cosmic origins
Astronomers thought they’d finally figured out where gold and other heavy elements in the universe came from. In light of recent results, they’re not so sure.
Image credit: Science Source/Tom McHugh.

Similar Articles