Soil bacterial communities vary more by season than with over two decades of experimental warming in Arctic tussock tundra

Climate change is altering the Arctic with unprecedented speed; the Arctic is warming at a rate two to six times the global average (Cohen et al., 2014). This has led to permafrost thaw and a deepening active layer within Arctic soils (Yi et al., 2018; Plaza et al., 2019). Concurrently, Arctic tundra systems are experiencing an increase in both overall plant productivity and shrub dominance (Tape et al., 2006; Wang et al., 2019)—a “greening of the Arctic” (Jia et al., 2003). These coupled aboveground-belowground changes in response to warming are expected to significantly alter tundra biogeochemical cycling and C balance. However, it remains unclear how these observed and predicted ecosystem-level changes correlate to, or are driven by, changes in the soil bacterial community (Virkkala et al., 2019).

Warming leads to shrubby expansion in tussock tundra ecosystems (Myers-Smith et al., 2011), with varied influences on the soil environment. On the one hand, the greater shade provided by shrubs can reduce summer soil temperatures and reduce the thickness of the seasonally thawed active layer (Blok et al., 2010). On the other, shrubs can enhance summer warming by increasing evapotranspiration and lowering albedo (Bonfils et al., 2012). Shrubs can also increase winter soil temperatures by trapping insulating snow (Loranty and Goetz, 2012). Microbes remain active in the soil under subzero temperatures under the snowpack (McMahon et al., 2011); because microbial activity is highly sensitive to temperature changes in frozen soils (Q₁₀ values near 10), even slight warming can greatly accelerate microbial activity (Mikan et al., 2002), making winter biogeochemical processes sensitive to changes in temperature. For instance, increased winter soil temperatures enable higher rates of nutrient mineralization (Schimel et al., 2004), leading to greater inorganic nutrient availability in spring (Mörsdorf et al., 2019), but also greater carbon starvation of the microbial community by late winter (Buckeridge and Grogan, 2008). Microbial communities within shrub soils also appear less susceptible to freeze-thaw induced lysis than those in Eriophorum vaginatum tussocks (Sistla et al., 2019). Nonetheless, shrub bacterial communities show greater seasonal variation than do tussock tundra communities, where active microbes are indistinguishable between late summer and early winter (McMahon et al., 2011). Together, this indicates that microbial communities may change drastically over winter under climate warming.

In addition to directly stimulating microbial decomposer activity, climate warming may alter microbial community structure through changes in the quality of plant inputs. Betula nana shrubs produce more labile root exudates than Eriophorum vaginatum tussock-forming sedges they replace (Proctor and He, 2017; Lynch et al., 2018), which may alleviate labile C limitation of microbial communities (Lynch et al., 2018). Indeed, bacterial phyla such as Acidobacteria that are associated with low resource availability (Fierer et al., 2007) were more abundant in tussock soils than in shrub soils (Wallenstein et al., 2007; Shi et al., 2015). On the other hand, Betaproteobacteria associated with high resource availability were more abundant in shrub soils (Wallenstein et al., 2007; Shi et al., 2015). Bacterial taxa associated with more copious and labile root exudation may be characterized as copiotrophs (Fierer et al., 2007; Nemergut et al., 2016). Greater labile C availability can also stimulate increased competition for nitrogen (N) between plants and microbes (Månsson et al., 2009), which may explain the increased abundance and diversity of nitrogen fixers under experimental warming (Feng et al., 2019). However, the effect of warming on plant communities is just one dimension of how warming changes the soil environment, and microbial community warming response may not parallel that seen between sedge- and shrub-dominated sites.

Bacterial communities show mixed responses to experimental warming in tundra ecosystems. For instance, bacterial community evenness declined in both surface and subsurface soils following 18 years of summer warming in a moist acidic tundra ecosystem (Deslippe et al., 2012). However, 15 years of experimental summer warming in a tundra heath ecosystem only affected microbial community structure in the surface soil, while the Gram positive to Gram negative bacterial ratio was unaffected by warming in both surface and subsurface soils (Rinnan et al., 2007). Microbial communities may also shift in response to warming on seasonal time scales (Buckeridge et al., 2013). However, this seasonal response is not always apparent (Deslippe et al., 2012), possibly because of slow turnover of identifying markers which obscures changes in the active community (Carini et al., 2020). Therefore, approaches that capture both actively growing and total microbial communities are best suited to understand how they change in response to warming or other environmental perturbations on short- to intermediate time frames.

To characterize how 22 years of experimental summer greenhouse warming affected the active layer bacterial community, we sequenced the actively growing and total communities across five seasonal time points in an Arctic tussock tundra system. Deciduous shrubs have approximately doubled in biomass with the warming treatment at this site, which allows the accumulation of additional insulating snow and therefore higher over-winter soil temperatures (Sistla et al., 2013). Soils derived from greenhouse plots also showed greater temporal variation in their potential hydrolytic extracellular enzyme activity during the year these soils were sampled, with the greatest increase in activity due to greenhouse treatment occurring during thaw and in the mineral soil (Sistla and Schimel, 2013). Nutrients and microbial biomass are temporally dynamic and peaked during late winter/thaw but were not found to vary between greenhouse treatments. In the current study, we hypothesized that experimental greenhouse warming—and the concurrent shifts in vegetation structure and changes in quantity and quality of inputs to soil—would lead to stronger effects on bacterial community structure than seasonality. This led to the specific predictions that: (1) the effect of warming on bacterial communities would be particularly apparent in late winter, (2) the effect of warming would be more readily apparent in the actively growing compared to the total bacterial community, (3) bacterial communities would diverge most rapidly during the transition from frozen to unfrozen soils, and (4) temporal changes in bacterial community would be less pronounced in mineral compared to organic soil due to dampened temperature oscillations in the former.

Our focal soils were collected from a long-term experimental warming study at the Toolik Long Term Ecological Research Site in Alaska (68°38′N, 149°34′W). The soil is a coarse-loamy, mixed, acidic, gelic Typic Aquiturbel, with a 30–50 cm thick organic horizon underlain by silty mineral soil (Romanovsky et al., 2007). The field experiment was initiated in 1989 and is comprised of four blocks of spatially paired plots with a greenhouse warming treatment and a control, unmanipulated plot. Passive summer warming was created by erecting clear polyethylene-sheets over permanent wooden frames (2.5 × 5 × 1.5 m) when the ground is snow-free (approximately June 1). This passively increases air temperature inside the greenhouses by an average of 2.1 °C. Although the effect of the greenhouses on summer soil temperatures appears to have diminished through time—possibly due to shading from the taller shrubby vegetation—the effect on winter soil temperatures has increased in the mineral soil, likely due to increased trapping of snow by the shrubs (Sistla et al., 2013). Seasonal maximum thaw depth at the site varies from 35–63 cm and is on average 22% greater in the greenhouses (55 cm vs. 45 cm; Sistla et al., 2013). The greenhouse warming treatment did not detectably affect soil moisture (Deslippe and Simard, 2011; Sistla et al., 2013), and air circulation occurred beneath the greenhouse bases due to uneven microtopography (Clemmensen and Michelsen, 2006). The greenhouse reduced photosynthetically active radiation and direct precipitation inputs but does not negatively influence plant growth (Deslippe and Simard, 2011). Vegetation at the site was initially dominated by the tussock-forming sedge Eriophorum vaginatum, but Betula nana shrubs had come to dominate the warmed plots by the time of soil collection, increasing from 16% cover in 1998 to 31% in 2008 (Sistla et al., 2013).

Soils were collected from the plots at five times during 2010, corresponding to late winter (April 28), thaw (May 19), summer (July 3), senescence (September 10), and early winter (November 10), as described in Sistla and Schimel (2013). The soil was frozen during the April, May, and November samplings, and so soils were collected using a hammer and chisel. After removing visible leaf litter, the soil was separated into upper organic (0–5 cm), lower organic (>5 cm), and mineral soil (between 5 and 10 cm below the organic, depending on variation across plots). The organic horizon contains a dense mat of roots, and so the soil can be considered root-influenced (“rhizospheric”) throughout. The soil was shipped back to UC Santa Barbara. Attempts were made to keep soils at temperatures close to those at which they were collected, from sampling until the end of incubations (see below). Therefore, soil was shipped frozen when soil was frozen in the field (April, May, and November), and at 4 °C in July and September when soils in the field were thawed. Descriptions of potential α-glucosidase, β-glucosidase, β-xylosidase, N-acetylglucosaminidase, acid phosphatase, phenol oxidase, and peroxidase enzyme activities, and soil carbon, nitrogen, and microbial biomass data can be found in Sistla and Schimel (2013), where they were originally published.

To explore how actively growing bacterial communities responded to the greenhouse treatment across soil depth and time of sampling, we incubated 5 g field moist soil subsamples in 100 ml falcon tubes with bromodeoxyuridine (BrdU). BrdU is a thymidine analog that can be used to analyze the proliferation of cells due to its incorporation into newly synthesized DNA (Borneman, 1999; Goldfarb et al., 2011; McMahon et al., 2011; Evans et al., 2014). Soil manipulations were completed in a 4 °C walk-in to help maintain soil temperatures during processing, and assays were initiated within 48 h of arrival in the lab (within a week of soil collection given transit time from Alaska). Soils were hand homogenized by horizon at the block level, and live coarse roots and rocks were removed. BrdU was added dropwise to soil samples as 1 ml of a 1.7 mM BrdU solution (made in sterile, Mili-Q H2O); the no BrdU control had 1 ml water added instead. Immediately following the BrdU or water-control addition, soils were gently vortexed for 20 s, which allowed the solution to visually mix into the soil. Soils were then incubated at the approximate field temperature at the time of soil collection, for the duration of time shown to be sufficient for the BrdU label to be detected in the DNA of incubated samples based on a preliminary study from greenhouse soils collected in May 2009 and other Arctic tundra studies (McMahon et al., 2009, 2011). Soils collected in April were incubated for 2 weeks at –10 °C, soils collected in May and November were incubated at –4 °C for 2 weeks, and soils collected in July and September were incubated at 1 °C for 10 days. Since soils were incubated at temperatures representative of those experienced in the field at the time of collection, data represent the communities active at that temperature at the time of collection rather than shifts in the taxa that would be expected to be active under a common laboratory temperature. Full information on soil handling temperatures can be found in Table 1.

Immediately following the incubation, DNA was extracted using MoBio PowerSoil Kits (MoBio Laboratories, Carlsbad, CA, USA). Antibody immunocapture was completed on a subsample of total DNA to isolate labeled DNA (containing the BrdU molecule) from unlabeled DNA extracted from a sample incubated only with water (McMahon et al., 2011). One unlabeled sample was also run with each immunocapture run to verify that there was no off-target capture of DNA. To check that BrdU immunocapture was effective, we amplified the immunocaptured DNA from both labeled and unlabeled samples using PCR with the primers 8f (5′-AGAGTTTGATCCTGGCTCAG; Turner et al. 1999) and 1389r (5′-ACGGGCGGTGTGTACAAG; Osborn et al. 2000), and then ran the PCR products on an agarose gel. The immunocapture run was deemed successful and on-target if the BrdU labeled DNA generated a band on the gel and no band was visible for the unlabeled immunocaptured control. Extracted DNA was stored at –20 °C for approximately 3 years prior to sequencing.

The V4 region of the 16S rRNA gene was amplified in triplicate using the original Caporaso primers 515F (GTGCCAGCMGCCGCGGTAA) and 806R (GGACTACHVGGGTWTCTAAT; Caporaso et al., 2011). Some samples had insufficient DNA for triplicate amplifications, and low sequence yield excluded them from downstream analysis. PCR products were cleaned with a QIAquick™ PCR Purification Kit (Qiagen, Valencia, CA, USA), and subsequently subject to sequencing using a 2 × 150 bp MiSeq run. This generated 16,018,656 paired-end reads of 151 bp prior to filtering and demultiplexing.

Sequence data processing was completed in QIIME2 version 2019.10 (Bolyen et al., 2019). Demultiplexing led to a total of 6,683,981 reads (median: 31,714, range: 1–82,661 per sample). Merging of forwards and reverse reads was completed using VSEARCH (Rognes et al., 2016) with minimum and maximum quality parameters set to 0 and 41, respectively. Merged paired-end reads were quality filtered based on PHRED score to trim when greater than three consecutive nucleotides had a PHRED score 4 or lower, and to exclude reads with 1+ ambiguous nucleotides. We used deblur (Amir et al., 2017) to remove erroneous reads and sort the reads into sub-OTUs (error-checked unique sequences), using trimmed read length of 250nt, keeping only sub-OTUs with at least 10 reads across all samples and found in at least two samples. This resulted in a total of 1,500,984 reads across 204 samples, with a range of 3–18,622 reads per sample. Taxonomy was assigned using a classifier trained to the V4 region of the 16S rRNA gene of the greengenes 99% identity reference sequences, with a confidence cutoff of 0.7 (Pedregosa et al., 2011). A phylogenetic tree was built by inserting the reads into the greengenes v. 13.8 reference tree using SEPP, which has been shown to provide more accurate placement for short reads characteristic of amplicon data than de novo-based tree building methods for sub-OTUs (Janssen et al., 2018). Sub-OTUs that could not be placed in the tree were removed from analysis. Reads classified as mitochondria and chloroplasts were subsequently filtered out of the tree, and the OTU table was rarefied to 2,250 reads per sample. This led to a total of 197 samples and 7,606 sub-OTUs for downstream analysis.

A single outlier (sample 268) was removed from downstream analysis based on substantially lower alpha diversity metrics (P value of χ² statistic > 0.999). The dominant sequence in this sample was very distant from the closest sequence in the database (>10 nearest sequenced taxon index (NSTI)), and so was likely to be a chimera. A final tally of replicates for each time point can be found in Table 2.

Weighted UniFrac distance (Lozupone and Knight, 2005) was calculated using functions in QIIME2. Differentially abundant dominant phyla were identified using ANCOM with plot as a random effect and the criteria that at least 70% of the comparisons were different (W = 7) and a P-value threshold of 0.05 (Mandal et al., 2015). All other downstream analyses and statistical tests were completed using R v. 3.6.1 after converting the QIIME output files to R-readable files with the qiime2R package v0.99 (Bisanz, 2018; R Core Team, 2019). We were generally interested in how actively growing and total communities differed from one another and how communities differed between soil horizons, experimental warming, and timing of sampling. Experimental factors driving beta diversity were tested with PERMANOVA using the adonis function in the R package vegan v. 2.5–6 (Oksanen et al., 2019), where permutations were completed within experimental block to account for nonindependence of warming treatment within block, sampling time within plot, and horizon within plot. We calculated the marginal sum of squares here due to imbalances in our sampling design.

The betadisper function was used to confirm that changes in community structure between treatments were real, and not simply because within-group differences under one level of a factor were larger than within-group differences under another level of the factor (Anderson, 2006). Therefore, this analysis was only run when adonis indicated differences between groups. Continuous environmental factors associated with community structure were assessed using distance-based redundancy analysis (McArdle and Anderson, 2001; Geml et al., 2015). We removed variables with variance inflation factors greater than 5 from our dbRDA analysis, keeping the variable we had the most complete data set for and/or was the most integrative metric (e.g., total extracellular enzyme activity rather than the activity of individual enzymes). We also calculated the Spearman correlation coefficients between relative abundance of dominant phyla and environmental variables, using the Benjamini-Hochberg correction to adjust for multiple testing (Benjamini and Hochberg, 1995).

Sparse subOTU matrices prevented the use of taxon-based measures of community turnover, so we used a phylogenetic beta diversity metric and a metric of change in phylum abundance. For the first, we calculated pairwise UniFrac distance between the same horizon and plot for successive time points divided by the number of days between samplings. This is a metric of temporal beta diversity (Buckley et al., 2019; Magurran et al., 2019), equivalent to the pairwise dissimilarity between communities collected from the same plot at consecutive time points. The second method was phylum-level temporal beta diversity of communities, which we calculated as Hedges G (Hedges and Olkin, 2014) using the effsize package (Torchiano, 2020). This metric accounts for both mean change and variation in the change across samples, such that a phylum which increases in relative abundance between two time points in all plots will receive a large positive score, and those which consistently decrease will receive a large negative score. We did not standardize Hedges G by time since the scientific community is generally familiar with it in its original units. Significant differences in temporal beta diversity were calculated using linear mixed models (lmer function in the lme 4 package v. 1.1–21; Bates et al., 2015)) followed by examination of specific comparisons of interest using the multcomp package (v. 1.4–12; Hothorn et al., 2008)). When fitting mixed models, we nested soil sample within plot and used experimental block as a random effect. This meant that all comparisons between actively growing and total communities were completed using paired total-active communities from the same soil sample, and all horizon effects were determined by comparing soils collected from the same plot on the same date. All figures were generated using ggplot2 v. 3.2.1 (Wickham, 2016).

Across greenhouse treatment, sample period, and depth, total communities were dominated by bacteria; archaea were absent from most samples (mode: 0, median: 0, range: 0%–0.88% of reads). Acidobacteria and Proteobacteria (13% of class-level reads were Alphaproteobacteria, 6.8% Betaproteobacteria, 4.0% Deltaproteobacteria, and 7.9% Gammaproteobacteria) dominated the communities (Figure 1). ANCOM indicated that Chloroflexi and Gemmatimonadetes were most abundant in the deeper mineral soil (Figure 1), while Armatimonadetes was more dominant in the organic soil (W = 8 in all cases). ANCOM did not identify any of the dominant phyla to differ in abundance under warming (W = 0 in all cases).

Twenty-two years of passive greenhouse warming did not alter tundra active layer bacterial community structure, PERMANOVA pseudo-F(1, 111) = 0.835, P = 0.84; dispersion pseudo-F(1, 111) = 0.013, P = 0.91; Figure 2. Instead, bacterial community structure was most strongly driven by horizon, R² = 0.206, pseudo-F(2, 111) = 16.26, P = 0.001; dispersion pseudo-F(2, 111) = 12.65, P < 0.001; mineral > shallow organic = deep organic, and time of sampling, R² = 0.128, pseudo-F(4, 111) = 5.04, P = 0.001; dispersion pseudo-F = 7.30, P < 0.001, September > April = May = July = November, with the horizon effect apparent in May, R² = 0.248, pseudo-F(2, 20) = 2.96, P = 0.001, dispersion pseudo-F(2, 20) = 0.636, P = 0.54, July, R² = 0.410, F(2, 21) = 6.61, P = 0.001; dispersion pseudo-F(2, 21) = 2.61, P = 0.10, September, R² = 0.476, pseudo-F(2, 21) = 8.66, P = 0.001; dispersion pseudo-F(2, 21) = 1.25, P = 0.309, and November, R² = 0.322, pseudo-F(2, 22) = 4.74, P = 0.001; dispersion pseudo-F(2, 22) = 1.40, P = 0.269. The sampling time point effect was particularly evident in the mineral soil (R² = 0.396 vs. 0.272 and 0.226 for surface and deep organic, respectively; dispersion pseudo-F(4, 39) = 1.84, P = 0.14, pseudo-F(4, 34) = 1.12, P = 0.36, and pseudo-F(4, 36) = 0.81, P = 0.528. Greenhouse also plots did not show stronger temporal variation in their total composition than control plots (rate of change in UniFrac distance between samples collected from the same plot at consecutive time points did not differ between greenhouse and control plots for any timeframe by soil horizon; Figure 3). Therefore, we did not distinguish between greenhouse treatments to evaluate how communities changed between time points or correlated with environmental factors, below.

Total bacterial communities are inferred to have diverged from one another at the lowest rate over winter (based on the difference between November and April communities) and at the greatest rate between April and May (Figure 3). These temporal shifts in community structure were not associated with strong changes in the relative abundance of any particular dominant phylum. Indeed, the only consistent phylum-level patterns observed were that Acidobacteria in the mineral soil were least abundant in September (Figure 1), and Actinobacteria decreased in relative abundance in the organic soil between May and July (Figure S3). This indicates that some dominant phyla showed a horizon-specific temporal pattern. The phylogenetic temporal beta diversity rate (rate of change in UniFrac distance) was similar among horizons, except for between September and November, when it was greater in the mineral soil compared to the surface and deep organic soils (Figure 3).

Total bacterial community composition was most strongly correlated with total organic carbon (R² = 0.63) and microbial biomass (R² = 0.31 based on dbRDA). The relative abundances of Chloroflexi, –3.95 ± 0.38 SE, t(102) = –10.46, P < 0.001, and Gemmatimonadetes, –2.55 ± 0.24, t(102) = –10.62, P < 0.001, were negatively correlated with total organic carbon, while the relative abundances of Actinobacteria, 1.82 ± 0.57, t(102) = 3.19, P < 0.05, and Armatimonadetes, 0.47 ± 0.08, t(102) = 5.89, P < 0.001, showed a positive correlation (Figure S4). Chloroflexi and Gemmatimonadetes relative abundance were negatively correlated with potential hydrolytic extracellular enzyme activity, while Armatimonadetes relative abundance was positively correlated (Figure 4A). Overall bacterial community structure also loosely correlated with total enzyme activity (R² = 0.11, P < 0.05), but did not correlate with microbial biomass-standardized total hydrolytic extracellular enzyme activity (R² = 0.00, P = 0.93). The relative allocation to the six hydrolytic enzymes assayed was correlated with total microbial community structure (Procrustes ordination of UniFrac distance vs. Hellinger distance of extracellular enzyme activity, m12 squared = 0.8, correlation = 0.45, P < 0.001).

Actively growing bacterial communities were substantially enriched in Bacteroidetes compared to total communities overall (24% higher; Figure S1; ANCOM with paired IC data, Wilcoxon test P < 0.05). However, a difference in bacterial community structure between total and actively growing communities was only apparent in July and September, when the soils were incubated at above-zero temperatures; Figure S2; time of sampling × activity interaction R² = 0.138, pseudo-F(4, 194) = 12.7, P = 0.001. Samples differed in their subOTU composition for both actively growing and total communities; therefore, the actively growing community was not a subset of the total community. Changes in actively growing BrdU-labeled communities were driven primarily by sampling time point, R² = 0.541, pseudo-F(4, 82) = 26.21, P = 0.001; dispersion pseudo-F(4, 82) = 6.77, P < 0.001; September = July > May = November, and, to a much lesser degree, horizon, R² = 0.069, pseudo-F(2, 82) = 6.69, P = 0.001; dispersion pseudo-F(2, 82) = 0.028 P = 0.97. There was no effect of the greenhouse treatment on actively growing community structure, R² = 0.004, pseudo-F(1, 82) = 0.850, P = 0.47; dispersion pseudo-F(1, 82) = 0.368, P = 0.546. We note, however, that we had insufficient replication to reliably assess changes in the actively growing community in the April and November samplings (Table 2) due to low yield of labeled DNA and poor sequencing for the actively growing community. Rate of change in bacterial community composition was greater in actively growing compared to total communities for both time periods we had sufficient replicates to examine approximately 45% higher in between July and September, t(16) = 6.13, P < 0.001, and 2.32× as high between May and July, t(10) = 14.65, P < 0.001, across all soil horizons (Figure 3).

Based on redundancy analysis, actively growing community structure correlated more strongly with ammonium concentration (R² = 0.46, P = 0.001) than with microbial biomass or organic carbon and nitrogen availability (R² = 0.28 – 0.38, P < 0.01). The actively growing community was only weakly correlated with total extracellular enzyme activity (R² = 0.15 P < 0.01 and R² = 0.05 ns for total and mass-specific enzyme activity). Paralleling patterns in the total community, the relative abundance of Gemmatimonadetes was negatively correlated with all variables except potential peroxidase activity and soil ammonium concentration (Figure 4B). Proteobacteria and Planctomycetes were both present at greatest relative abundance in samples with high TOC and low potential peroxidase activity, and in soils high extractable nutrient concentrations and potential hydrolytic extracellular enzyme activity (Figure 4B, S5). Verrucomicrobia were also present in greater relative abundance in the actively dividing communities when ammonium and extractable organic carbon were high, but their abundance did not correlate with extracellular enzyme activity (Figure 4B). The relative abundance of Chloroflexi and Gemmatimonadetes in the actively dividing community data decreased with TOC concentration (Figure S5) and with the potential activity of β glucosidase, N-acetylglucosaminidase, and acid phosphatase enzymes (Figure 4B).

Numerous studies have found that warming alters Arctic plant community structure and productivity, typically increasing woody plant dominance. However, how warming alters tundra soil microbial communities and what role vegetation changes play in driving bacterial community shifts have remained less clear. We found bacterial community structure was unaffected by 22 years of passive summer warming despite greenhouse treatment increasing woody plant dominance. However, both total and actively growing soil bacterial communities were temporally dynamic, indicating that they may be more sensitive to short-term seasonal changes than long-term changes in mean temperature.

We hypothesized that bacterial communities turn over most rapidly during the transition from frozen to unfrozen soils (May–July), which was observed for the actively growing, but not total, communities. Microbial community composition has been observed to be temporally dynamic at other tundra sites (Wallenstein et al., 2007; McMahon et al., 2011; Vořšková et al., 2019), and—in contrast to a previous study (Deslippe et al., 2012)—our results indicate that bacterial communities from this long-term warming study are also temporally dynamic. Rapid microbial community change may occur predominantly when soil thaws in spring because the labile C, N, and phosphorus (P) resulting from freeze-lysis of microbes in fall becomes more available for the surviving microbes to take up (Schimel and Clein, 1996; McLaren et al., 2017; Sistla et al., 2019). Concurrently, plant roots and their mycorrhizae track the nutrient-rich thaw front in soil (Hewitt et al., 2019, 2020), which also presumably provides fresh labile C inputs to the surrounding microbial communities in the deeper soil (Iversen et al., 2015). Therefore, it is not surprising the rate of change in community structure in the actively growing community was greatest between May and July. Total bacterial communities shifted most rapidly between April and May, however, which may be due to increased liquid water availability as the active layer begins to thaw (Ishizaki, 1994; Wen et al., 2012). Curiously, the high community temporal beta diversity observed between April and May was not apparent at the phylum level for the total communities. This indicates that despite strong correlations between soil variables and relative abundance for some phyla, OTUs within a phylum do not act in concert and/or the sampling time point effect was heterogeneous across plots and horizons.

Given known influences of chronic warming on plant community and therefore the potential for the soil to maintain a warm deeper insulating layer of snow over winter, we also anticipated that bacterial communities would be more temporally variable in control plots. However, we did not find an effect of warming on temporal beta diversity rate, contrasting with both results showing greater bacterial community temporal beta diversity in shrub compared to tussock soils (McMahon et al., 2011), and increased bulk community temporal beta diversity under freeze/thaw in Eriophorum vaginatum soils (Sistla et al., 2019). These findings suggest that changes in microbial communities may not always follow changes in aboveground plant communities—as has been previously reported for bulk soil (Singh et al., 2007)—as warming acts directly on microbes in addition to indirectly through plant inputs.

In contrast to some other soil warming studies (Feng et al., 2019; Johnston et al., 2019; Vořšková et al., 2019), we did not find 22 years of experimental warming had changed bacterial communities. This is despite evidence for warming-driven changes in the microbial eukaryote community (Sistla et al., 2013) and extracellular enzyme activity (Sistla and Schimel, 2013) at this site. Shrub-dominated soils have been observed to be enriched in Proteobacteria and depleted in Acidobacteria compared to Eriophorum-dominated ones (Wallenstein et al., 2007), so the mechanisms by which warming-driven shrub expansion occurred but the tussock bacterial community persisted remains to be explored. One possibility is that changes in plant communities are primarily associated with shifts in fungal rather than bacterial communities, as mycorrhizal fungi provide up to 85% of the N to plants in these systems (Hobbie and Hobbie, 2006). Alternatively, rapid succession of plant communities in-situ may not drive changes in microbial communities in the same way that preexisting differences in abiotic conditions that allow for preferential establishment of Eriophorum sedges or Betula shrubs do. The soils and microtopography of greenhouse plots have largely retained the physical structure of tussock tundra at our site (Wallenstein et al., 2007; Sistla et al., 2019), and dissolved C and N availability were largely unaffected by warming in soils collected for this study (Sistla and Schimel, 2013). Nonetheless, based on automated ribosomal intergenic spacer analysis, bacterial and fungal communities had shifted over 18 years of warming during growing season sampling at this tundra site (Deslippe et al., 2012). The divergence between these previously reported results and our own may reflect differences in the resolution and taxonomic bias of methods used (Wallenstein et al., 2007; Campbell et al., 2010; Kovacs et al., 2010; Kim et al., 2014; Ricketts et al., 2016) or because the effect of warming on microbial communities is inconsistent through time (Melillo et al., 2017). Finally, it is possible that storing and incubating soils from greenhouse and control plots at the same temperature and without plant inputs for 2+ weeks may have annulled any impacts warming may have had on microbial communities in the field, which could have been observed if DNA had been extracted from field-fresh soils.

Compared to previous studies of microbial communities in this warming experiment, our phylotype-level data enabled us to evaluate the phylogenetic structure of the bacterial community and which bacteria dominate moist acidic tussock soils through the active layer and across time points. Further, by comparing this community data to concurrently measured edaphic factors, we were able to determine which soil variables most strongly correlate with overall bacterial community structure. For instance, Chloroflexi and Gemmatimonadetes, physiologically diverse groups previously observed to increase in dominance with depth (Kim et al., 2016), correlated positively with soil oxidative enzyme activity but negatively with hydrolytic enzyme activity. A positive association between the relative abundance of these two phyla and lignin breakdown has been previously been observed in aqueous environments (Colatriano et al., 2018; Zhang et al., 2019) and may indicate that these taxa favor the consumption of lignin decomposition products over polysaccharides in these soils.

Procrustes analysis confirmed that the bacterial phylogenetic community structure and whole soil community potential extracellular enzyme activity matrices were correlated, indicating that different bacterial communities may preferentially target certain substrates for decomposition. However, the overall extracellular enzyme activity was largely a function of resource availability and microbial biomass, as total hydrolytic extracellular enzyme activity correlated with the community ordination axes, but mass-specific extracellular enzyme activity did not. This indicates that bacteria in high and low biomass communities allocate similar effort overall to hydrolytic extracellular enzyme activity but differ in how that effort is allocated to different enzymes. This does not preclude the possibility that it is actually fungi, rather than bacteria, who are primarily important for producing extracellular enzymes in these soils (Schneider et al., 2010), such that the relationship between enzyme classes and bacterial community is purely correlative.

Microbes in high latitude tundra ecosystems are thought to be highly susceptible to climatic change due to the large degree of environmental filtering that occurs during community assembly (Bahram et al., 2018). Although two decades of experimental summer warming shifted the plant community from tussock to shrub dominance, this warming did not alter bacterial community structure in the active layer. However, both the total and actively growing communities differed between time points which differed in their temperature and plant phenological phase. This pattern indicates that these communities were more directly responsive to changes in temperature than the long-term changes in the soil environment such as in the structure of plant and protist communities (Sistla et al., 2013). Furthermore, the stability of the microbial community in response to over two decades of warming suggests that changes in large-scale biogeochemically relevant processes such as net primary productivity, microbial biomass carbon, and respiration rate may occur even without substantial changes in the soil bacterial community structure.

The sequencing data used in this study are available in the SRA under accession numbers SAMN14849733-9959. The scripts used to create the figures used in this analysis can be found at https://osf.io/kav4f/.

The supplemental files for this article can be found as follows:

Figure S1. Phylum-level relative abundance of reads in immunocaptured and total 16S communities. Only samples for which paired total-immunocaptured data were available are included (n = 75 pairs). ANCOM was used to determine which phyla differ between immunocaptured and total communities, using a threshold of 0.8 and a P-value cutoff of 0.05.

Figure S2. Unifrac-based principal coordinates analysis plot of bacterial communities. A single ordination was generated for all samples, but samples are separated out into different panels according to time of sampling and soil horizon. Symbol color denotes whether the sample was derived from a greenhouse or control plot, and symbol shape denotes whether the sample came from a BrdU ("actively growing") or the total community.

Figure S3. Change in relative abundance of dominant phyla between timepoints. Data are presented with one line per phylum, separated by soil horizon sampled. Only the total communities are shown. The metric depicted is Cohen's D, calculated as the difference in mean abundance between the two time points, divided by the pooled standard deviation of relative abundance weighted by the sample size. Thus, the value is an uncertainty-standardized measure of change in relative abundance for the phylum. Rates are not corrected for the time between successive samplings.

Figure S4. Relationship between the concentration of total organic carbon (here shown as the log) and relative abundance of dominant phyla in the total communities. Plot was included as a random effect in a linear mixed model to account for sampling horizons within plots and plots multiple times across time points. Only regressions where the Benjamini-Hochberg-corrected P-value for the slope being different from zero is less than 0.05 are shown.

Figure S5. Relationship between the concentration of total organic carbon (here shown as the log) and relative abundance of dominant phyla in the actively growing communities. Plot was included as a random effect in a linear mixed model to account for sampling horizons within plots and plots multiple times across timepoints. Only regressions where the Benjamini-Hochberg-corrected P-value for the slope being different from zero is less than 0.05 are shown.

The authors wish to thank Noah Fierer’s lab for completing the sequencing used in this study.

Funding for this project came from a National Science Foundation Doctoral Dissertation Improvement Grant (Award # 1110843) and an NSF EAGER grant to S.A.S. (Award # 1841610).

The authors declare no competing interests related to this manuscript.

Contributed to conception and design: SAS, JPS.

Contributed to acquisition of data: SAS.

Contributed to analysis and interpretation of data: GP.

Drafted and/or revised the article: GP, SAS, JPS.

Approved the submitted version for publication: GP, SAS, JPS.

How to cite this article: Pold, G, Schimel, JP, Sistla, SA. 2021. Soil bacterial communities vary more by season than with over two decades of experimental warming in Arctic tussock tundra. Elementa Science of the Anthropocene 9(1) DOI: https://doi.org/10.1525/elementa.2021.00116.

Domain Editor-in-Chief: Steven Allison, University of California, Irvine, CA, USA

Associate Editor: Stephanie Kivlin, Ecology and Evolutionary Biology, University of Tennessee, Knoxville, TN, USA

Knowledge Domain: Ecology and Earth Systems