Impact of long-term fertilizer and summer warming treatments on bulk soil and birch rhizosphere microbial communities in mesic arctic tundra

ABSTRACT

Recent climate warming in the Arctic is enhancing microbial decomposition of soil organic matter, which may result in globally significant greenhouse gas releases to the atmosphere. To better predict future impacts, bacterial and fungal community structures in both the bulk soil and the rhizosphere of Arctic birch, Betula glandulosa, were determined in control, greenhouse summer warming, and annual factorial nitrogen (N) and phosphate (P) addition treatments twelve years after their establishment. DNA sequence analyses at multiple taxonomic levels consistently indicated substantial bulk soil and rhizosphere microbial community differences among the fertilization treatments but no significant greenhouse effects. These results suggest that climate warming will likely increase the activity rates of soil microbial decomposers but without substantially altering the structure of either the bacterial or fungal communities. Differential abundance testing revealed changes in ectomycorrhizal fungal species of the genus Thelephora in both bulk soil and rhizosphere, with increases in their relative abundance in P and N + P amended plots compared with warming and controls. Because birch is the principal low Arctic ectomycorrhizal host, our results suggest that these fungi may promote this shrub’s competitiveness where tundra soil nutrient availability is enhanced by warming or other means, ultimately contributing to arctic vegetation “greening.”

KEYWORDS:

• Climate change
• nitrogen and phosphate fertilization
• Betula glandulosa
• ectomycorrhizal fungi
• soil community diversity

Introduction

Arctic tundra ecosystems are characterized by their low productivity, associated with low nutrient availability, as well as low temperatures and a short growing season. The short summer thaw limits soil microbial activity and the associated supply of macronutrients, including nitrogen (N) and phosphorus (P), into plant available pools, which then restricts ecosystem productivity (Blaud et al. 2015). Arctic soils are also well known for their particularly high levels of accumulated soil carbon (C), which are a consequence of temperature-constrained microbial decomposition of organic matter. However, with increasing temperatures in the Arctic, microbial activity and rates of decomposition are expected to rise in tundra ecosystems (Shaver et al. 2006; Brzostek et al. 2012; Cavicchioli, Ripple, and Timmis et al. 2019). This, as well as increased depth of thaw into the permafrost, could increase overall organic matter decomposition and, as a consequence, the CO₂ release from the soil to the atmosphere will further contribute to climate warming (Schuur et al. 2009, 2013; Schaefer et al. 2014). Furthermore, the recent vegetation “greening” across much of the Arctic due to enhanced growth of deciduous shrubs has been specifically attributed to climate warming–induced increases in soil fertility (Chapin et al. 1995; Tremblay, Lévesque, and Boudreau 2012; Myers-Smith et al. 2015). Thus, there is high interest in better understanding the structure and functioning of arctic soil microbial communities and their responses to warming and increased nutrient availability.

Compared with other terrestrial ecosystems, the structure of bacterial and fungal communities in arctic soils has revealed a higher level of diversity than initially expected (Neufeld and Mohn 2005; Malard and Pearce 2018). It is well known that soil community microbiota can be influenced by vegetation type, as evidenced by comparisons of bacterial and fungal communities in arctic tussock, intertussock, and shrub sites (Wallenstein, McMahon, and Schimel 2007; Kumar et al. 2017). Similarly, bacterial communities differ among tundra wet sedge, birch hummock, tall birch, and dry heath environments, where distinct soil communities are also correlated with pH, ammonium concentration, moisture, and the C:N ratio (Martiny et al. 2011; Y. Shi et al. 2015). Although these correlations suggest that enhanced nutrient availability could alter microbial communities in arctic soils, the impacts of experimental fertilization have been more frequently investigated and are better understood in agricultural and pasture soil communities than in arctic soils (e.g., Jangid et al. 2008; Zhong et al. 2010; but see Deslippe, Egger, and Henry 2005; Prager et al. 2017).

Changes in tundra soil fertility resulting from increased seasonal temperatures are likely to be slow and cumulative (Ballhausen et al. 2020; Hobbie, Nadelhoffer, and Högberg 2002), which points to the importance of extended experiments to evaluate the effect of warming on soil communities. Long-term greenhouse experiments, as well as nutrient addition plots, were established in moist acidic tundra soils near Toolik Lake, Alaska, in 1988 as part of the Arctic Long-Term Ecological Research experiments to simulate aspects of prolonged climate change. After eighteen years, evaluation of the microbial communities in the warming plots indicated lower bacterial community evenness and greater fungal evenness in warmed than control plots (Deslippe et al. 2012). More specifically, the data suggested that the dominance of Actinobacteria was increased and Proteobacteria were reduced in the bacterial communities. Furthermore, the warming regime increased fungal abundance of taxa such as Cortinarius spp. and the ectomycorrhizal fungus Russula ssp. The effects of nutrient additions on microbial communities revealed that overall long-term combined N and P (N + P) fertilization resulted in significant changes in bulk soil microbiota, linked in turn to changes in C and N availability and aboveground plant communities (Campbell et al. 2010).

Another long-term summer warming and nutrient addition experiment was established at the Tundra Ecosystem Research Station at Daring Lake in the Northwest Territories in 2004, but in this instance for a mesic birch hummock ecosystem (Zamin and Grogan 2012). Here, greenhouse warming resulted in increased above- and belowground C, N, and P pools within the vegetation but did not increase microbial or soil solution C, N, and P pools. In contrast, after experimental fertilization, N and P amendments increased microbial pools of these nutrients (Zamin 2013). At that time, however, the impact of these experimental treatments on soil bacterial and fungal community structure was not examined. Furthermore, the microbial rhizosphere communities directly associated with or on the surfaces of plant roots was not explored, even though the rhizosphere is known to play distinctive and important roles in plant nutrient acquisition and health (Turner, James, and Poole 2013; Guyonnet et al. 2018). Here we separately characterized the microbial communities of the bulk soil and the rhizosphere of Arctic birch, Betula glandulosa, a predominant species in the low Arctic mesic tundra that appears to be expanding, possibly as a result of enhanced nutrient levels due to warmer soils and linked to climate change (Chapin et al. 1995; Zamin and Grogan 2012; Zamin, Bret-Harte, and Grogan 2014; Myers-Smith et al. 2015). Thus, our objectives were to evaluate the effects of long-term greenhouse summer warming and nutrient addition treatments on bacterial and fungal communities in the bulk soil and rhizosphere of B. glandulosa. We focused on B. glandulosa because it is a conspicuous sentinel of climate change in the Arctic and tested the hypothesis that twelve years of either summer warming or factorial N + P nutrient additions will significantly shift the bacterial and fungal communities of the bulk soil and Arctic birch rhizosphere.

Methods

Study site and experimental treatments

The experimental treatment plots were established in the summer of 2004 in mesic birch hummock tundra in the central Canadian low Arctic at the Tundra Ecosystem Research Station at Daring Lake, Northwest Territories (64°52′ N, 111°34′ W; Figure S1). Details of the treatment plots have been previously described (Zamin and Grogan 2012). Briefly, experimental plots (5 m × 7 m) were randomly allocated to one of the treatments, with five replicate plots per treatment (n = 5; Table 1). Control treatment plots were given no nutrient supplementation, nitrogen treatment plots were amended with 10 g N m⁻² yr⁻¹, phosphate treatment plots were modified with 5 g P m⁻² yr⁻¹, and N + P treatment plots were amended with both 10 g N m⁻² yr⁻¹ and 5 g P m⁻² yr⁻¹. N and P were applied in the form of regular agricultural-grade dry powdered NH₄NO₃ and triple superphosphate (45 percent P₂O₅), respectively. The summer warming experimental treatment was established using A-frame greenhouses covered with a heavy polyethylene film (0.15 mm) that was erected in early-mid-June and removed during late August to early September each year and that reduced mean daytime photosynthetically active radiation by 32 percent (Zamin, Bret-Harte, and Grogan 2014). Small vent holes (20 cm in height) in the apex of each end of the greenhouses were designed to reduce the possibility of extreme elevated temperatures and high humidity. Mean summertime air and upper soil (0–10 cm) temperatures averaged 2.5°C higher in the greenhouses (range 0°C–13.8°C) compared to adjacent control plots without greenhouses (Zamin and Grogan 2012). Multiple soil volumetric water content readings per plot on several summer sampling dates over the years of the study consistently indicated no significant effects of the greenhouse plastic covering on soil moisture (Zamin, Bret-Harte, and Grogan 2014), presumably because of subsurface lateral water flow from around the greenhouse edges. Nutrient analysis has been previously described and presented (Zamin 2013; Zamin, Bret-Harte, and Grogan 2014). It is important to note that all plots were initially chosen on the basis of their similarity in vegetation cover, soil type, and topography. Only then were they randomly allocated to the treatments and their respective controls, and thus our interpretation of the results necessarily assumes that the microbial communities in all plots were similar prior to the initiation of the experimental manipulations.

Table 1. Summary of annual experimental treatments in birch mesic tundra vegetation near Daring Lake that were initiated in 2004 in randomly allocated plots (5 m × 7 m) at five replicate plots per treatment and assayed in triplicate for this study

Treatment	Addition	Form added
Control	No addition	n/a
High N	10 g N m^−2 yr^−1	NH4NO3
High P	5 g P m^−2 yr^−1	TSP (45% P2O5)
High N + P	10 g N m^−2 yr^−1 and 5 g P m^−2 yr^−1	NH4NO3 and TSP
Warming	Polyethylene film installed from June to August	A-frame greenhouse

Sample collection and processing

Bulk soil and birch rhizosphere samples were collected on 30 August 2015 in triplicate from each of the five replicate plots within each experimental treatment (control, P, N, N + P, and warming; Table 1), to a total of n = 25. Three randomly selected birch (B. glandulosa) plants that were at least 1 m apart were selected in each treatment replicate plot. For the bulk soil samples, we first removed the overlying surface green moss tissue and leaf litter and then cut out an approximately 30 g fresh weight sample of organic soil (including brown moss) to approximately 5 cm depth from at least 25 cm away from the base of each one of the selected birch shrubs (n = 3 per plot). By contrast, the three birch rhizosphere samples from each plot were taken by uncovering and cutting out approximately 10-cm-long sections of root (totaling ~2 g fresh weight) that were clearly connected to the same birch plants described above. Both bulk soil and root samples were cut out from the ground using an ethanol-sterilized knife for each sample and transferred to sterile sample bags (Figure S1). All samples were kept on ice in an insulated cooler and transported to Queen’s University within five days, where they were then processed within 24 hours of arrival. In the laboratory, bulk soil samples were transferred to 15 mL conical tubes, stored at −20°C, and lyophilized prior to DNA extraction. For the rhizosphere samples, roots (~6.5 g) were taken from each sample bag and obvious soil particles adhering to the root sections were gently removed using 70 percent ethanol-sterilized scoops and forceps. The roots were then transferred into sterile 15 mL conical tubes and stored at 4°C until they were washed the following day.

Root samples were washed with 5 mL of cold 10 mM NaCl solution to remove soil and microbial cells in the vicinity of the root surface (Shakya et al. 2013). Tubes were then shaken by hand for 5 minutes and subsequently placed on a platform shaker at 4°C for an additional 30 minutes of shaking at 180 rpm. The roots were then carefully removed from each tube using sterile forceps, and the remaining saline solution and soil debris that had washed off the roots was stored at −20°C. To obtain the rhizosphere sample, the frozen saline wash solutions were thawed and centrifuged for 20 minutes at 4°C at 1500 × g (Bulgarelli et al. 2012). The supernatant was removed carefully to avoid disturbing the loose pellet. Centrifugation was repeated once more to remove as much supernatant as possible. The pellet was defined as the rhizosphere sample and stored at −80°C until DNA extraction.

DNA extraction

DNA was extracted from the lyophilized bulk soil samples (150 mg each) using the NucleoSpin Soil Extraction Kit (Macherey-Nagel; D-Mark Biosciences, Toronto, ON, Canada). The manufacturer’s protocol was followed with the addition of a repeated sample lysis step. Purified DNA was eluted with 50 μL of sterile double-distilled water (ddH₂O). Rhizosphere samples were also extracted using the NucleoSpin Soil Extraction Kit (Macherey-Nagel), following the recommended protocol but again with an additional sample lysis step, and the DNA was eluted with 30 μL of sterile ddH₂O. Sample quality and concentration were assessed using both agarose gel electrophoresis and a Nano-Drop 2000 spectrophotometer (Thermo-Fisher, Mississauga, ON, Canada).

DNA amplification and gene sequencing

All bulk soil and birch rhizosphere DNA extractions were sent to the Center for the Analysis of Genome Evolution & Function at the University of Toronto for sequencing of the V4 region of the 16S rRNA gene as well as a portion of the fungal internal transcribed spacer (ITS) region. The V4 hypervariable region was amplified using a universal forward sequencing primer (515 F) and a uniquely barcoded reverse sequencing primer (806 R) to allow for multiplexing (Caporaso et al. 2012). Amplification reactions for the soil samples were performed using 12.5 μL of KAPA2G Robust HotStart ReadyMix (KAPA Biosystems, Millipore-Sigma, Oakville, ON, Canada), 1.5 μL of 10 μM forward and reverse primers, 8.5 μL of sterile ddH₂O, and 2 μL of DNA. Amplification was achieved by cycling the reaction at 95°C for 3 minutes, 18 cycles of 95°C for 15 seconds, 50°C for 15 seconds, and 72°C for 15 seconds, followed by a 5-minute 72°C extension. All amplification reactions were done in triplicate alongside negative controls for each barcode, verified on a 1 percent agarose Tris-borate EDTA gel, and then pooled to reduce amplification bias. Pooled triplicates were combined at approximately equal concentrations. The final library was purified using 0.8 X magnetic Ampure XP beads, selecting for the bacterial V4 amplified band. The purified library was quantified and loaded onto the Illumina MiSeq for sequencing, according to manufacturer’s instructions (Illumina, San Diego, CA). Sequencing was performed using the MiSeq Reagent Kit v2 (150 bp X 2).

For fungal analysis, the ITS1 region of the 18S-ITS1-5.8S-ITS2-28S rRNA sequence was amplified using a ITS1F and ITS1R primer set (Willger et al. 2014). Amplification reactions were performed using 12.5 µL of KAPA2G Robust HotStart ReadyMix (KAPA Biosystems), 1.5 μL of 10 µM forward and reverse primers, 8.5 µL of sterile ddH₂O, and 1 µL of DNA. The ITS1 region was amplified by cycling the reaction at 95°C for 3 minutes, 25 cycles of 95°C for 15 seconds, 56°C for 15 seconds, and 72°C for 15 seconds, followed by a 5-minute 72°C extension. All amplification reactions were done in triplicate, verified on a 1 percent agarose Tris-borate EDTA gel, and then pooled to reduce amplification bias. Pooled triplicates were quantified using a Qubit fluorometer (Thermo-Fisher) and the standard Nextera XT protocol was followed to ligate adapters and 300 to 500 bp fragments were selected. The sample libraries were then quantified using fluorescence with a Qubit fluorometer and pooled. The final pooled library was quantified and loaded on to the Illumina MiSeq for sequencing according to manufacturer instructions (Illumina). Sequencing was performed using the MiSeq Reagent Kit v2 as previously noted.

Analysis of bacterial and fungal microbiota

For bacterial sequencing results, the UNOISE pipeline, available through USEARCH v9.2, was used for sequence analysis, and fungal sequencing was processed using the UNOISE pipeline, available through USEARCH v10.0.240 (Edgar 2010, 2013, 2016). The last base, typically error prone, was removed from all sequences. Sequences were assembled and quality trimmed using fastq_mergepairs and fastq_filter, with fastq_maxee set at 1.0 and 0.5, respectively.

For bacterial analysis, sequences less than 233 bp (20 bp shorter than the average length) were filtered from the data. Following the UNOISE pipeline, merged pairs were then de-replicated and sorted to remove singletons. Sequences were denoised and chimeras were removed using the unoise2 command. Assembled sequences were then mapped back to the chimera-free denoised sequences using 97 percent identity operational taxonomic units (OTUs). Taxonomy assignment was executed using utax, available through USEARCH, and the UNOISE compatible Ribosomal Database Project database v16, with a minimum confidence cutoff of 0.9 (Q. Wang et al. 2007). OTU sequences were aligned using PyNast accessed through QIIME (Caporaso et al. 2010). Sequences that did not align were removed from the data set and a phylogenetic tree of the filtered aligned sequence data was made using FastTree (Price, Dehal, and Arkin 2009).

For fungal ITS analysis, following the UNOISE pipeline, unique sequences were identified from the merged pairs. Sequences were denoised and chimeras were removed using the unoise3 command in USEARCH. Assembled sequences were then mapped back to the chimera-free denoised sequences at 97 percent identity OTUs using the otutab command. Taxonomy assignment was executed using SINTAX (Edgar 2016), available through USEARCH, and the SINTAX compatible Ribosomal Database Project Warcup ITS v2 database, with the default minimum confidence cutoff of 0.8 (Q. Wang et al. 2007).

Analysis of the generated OTU tables was done using QIIME v1.9.1 (Caporaso et al. 2010), including data normalization, taxonomic relative abundance, rarefaction, alpha diversity calculation, and beta diversity calculations. Data were normalized using both rarefaction to the lowest sample count and as well cumulative sum scaling (CSS) normalization (Paulson et al. 2013). Alpha diversity was calculated with both Faith’s phylogenetic diversity (PD) metric for bacterial communities and the Shannon index (also known as the H or Shannon-Wiener index or Shannon-Weaver index), using rarefied OTU tables. Beta diversity was calculated using the weighted UniFrac metric, a distance matrix, which takes into account the taxonomic relationship (i.e., phylogenetic closeness) as well as the abundance of OTUs within samples (Lozupone and Knight 2005). The distance matrix of beta diversity values for the samples was used as the input for nonmetric multidimensional scaling (NMDS) analysis based on Bray-Curtis to quantify compositional dissimilarity. and this was implemented through the vegan package (2.5–6) in R v3.3.2, with environmental variables (pH, NO₃-N, NH₄-N, and PO₄-P) collected from the same experimental plots (Zamin 2013). This then produced ordinations that could be displayed using the using the envfit function (vegan). As noted, calculations were performed on CSS normalization, based on the recommendations of McMurdie and Holmes (2014).

Statistical analyses

Differential abundance of bacterial OTUs at the family and phylum levels and fungal OTUs at the phylum and genus levels was tested using analysis of composition of microorganisms (ANCOM), implemented through QIIME2, using parametric one-way analysis of variance for hypothesis testing (Mandal et al. 2015). Differences in alpha diversity were tested using nonparametric analyses of variance, implemented through Prism 7. Analysis of the statistical significance of sample groupings for beta diversity metrics was done using the compare_categories.py script in QIIME, which uses the R vegan package to run both ANOSIM and adonis tests. Finally, statistical analysis of beta diversity distances within and between treatments was done using the make_distance_boxplots.py script in QIIME, with parametric t tests run including Bonferroni p value correction for multiple comparisons. Mantel correlation tests were performed using the script compare_distance_matricies.py in QIIME to compare beta diversity distances to soil chemical data including pH, NO₃-N, NH₄-N, and PO₄-P concentrations of samples collected from the same experimental plots in summer 2012 (Zamin 2013), with the assumption made that the relative differences in these values among plots did not change significantly in the intervening period. Prior to these tests, soil geochemical data were converted into a distance matrix using the script distance_matrix_from_mapping.py. All statistical analyses were interpreted using a significance level of p < .05.

One sample from the fungal rhizosphere sequencing data within the N + P treatment appeared to be an outlier and did not cluster within its own treatment or within any other treatment when visualized using NMDS. Thus, for analysis of both alpha and beta diversity, the sample was removed. However, it was retained for the testing of differential taxon abundance because there was no visible change in the average fungal relative abundance in the rhizosphere of the N + P treatment before and after removal of the sample.

Data deposition

Sequence data have been submitted to the National Center for Biotechnology Information with links to BioProject accession number PRJNA670726 in the NCBI BioProject database (https://www.ncbi.nlm.nih.gov/bioproject/).

Results

Data set reads

Tables generated from the V4 region of the 16S rRNA gene and ITS1 region sequences for the bacterial and fungal microbiomes, respectively, were produced using two different pipelines for both the bulk soil and birch rhizosphere sample groups, as detailed above. The bacterial bulk soil OTU table contained 874 unique OTUs in a total of 1,496,650 reads, averaging 59,866 reads per sample with a standard deviation of 35,608. The birch rhizosphere bacterial OTU table had relatively more unique OTUs at 1,679 in a total of 1,213,872 reads, with an average of 48,555 reads/sample (SD = 26,409). For the fungal sequencing, the bulk soil OTU table had 1,379 unique reads in a total of 787,869, thus averaging 31,515 reads/sample (SD = 17,135). The birch rhizosphere fungal sequences had fewer unique OTUs and total reads at 926 and 257,299, with an average 10,292 reads/sample (SD = 4,721).

Phylum-level effects of nutrient addition and warming treatments

The average relative abundances of both the bacterial and fungal phyla in the bulk soil and birch rhizosphere samples suggested differences between the experimental treatments and their respective controls (Figure 1). Most of such differences in the bacterial communities were observed in the N + P treatment plots compared with the control, unamended plots (Figure 1a). For example, there appeared to be a lower average abundance of the candidate division phyla WPS-1 and WPS-2 in the bulk soil of both the P and N + P treatments compared with controls, as well as a relatively lower average abundance of Acidobacteria in the N addition treatment plots (Figure 1a). Similarly, in the rhizosphere samples from the N + P addition, treatments there appeared to be a relatively lower average abundance in the phyla WPS-1 and WPS-2 compared with the control and warming plots (Figure 1a). Furthermore, the data suggested that Verrucomicrobia (Verrucomicrobiota) had a greater average relative abundance in the P addition treatment of bulk soil samples, with a lower abundance of Planctomycetes in the N + P treatment compared with the control and warming plots in both the bulk soil and birch rhizosphere (Figure 1a). Notably, there were minimal apparent differences in average relative abundances between the warming and control plots of bacterial phyla present in either the bulk soil or birch rhizosphere samples. To evaluate significant changes in the differential abundance of phyla, ANCOM was implemented. This analysis revealed that there was a significant difference in the differential abundance of the candidate phyla WPS-1 in both the bulk soil and rhizosphere, mirroring the changes in relative abudance observed between treatments (Table S1). Additionally, a significant differential abundance of the phylum Parcubacteria was detected in the rhizosphere, although this phylum had a low relative abundance (<1 percent; Table S1).

Figure 1. Average relative abundances of phyla in bulk soil and birch rhizosphere within each experimental treatment group for both the (a) bacterial and (b) fungal communities of low Arctic mesic birch hummock tundra. Sequence results analyzed to phylum level for (a) bacterial communities and (b) fungal communities in bulk soil and rhizosphere. The relative abundances shown represent the average of treatment replicates (n = 5), with treatment groups indicated on the x-axis. Any phyla present at less than 0.1 percent relative abundance across all treatment groups were pooled together.

Apparent visible differences in the average relative abundance of fungal phyla among treatments were also noted (Figure 1b). For example, the relative abundance of bulk soil and rhizosphere Basidiomycota appeared to be greater in both the P and N + P addition treatments compared with the control and warming plots (Figure 1b). Furthermore, the data indicated a decrease in the average relative abundance of taxa in the bulk soil assigned to “other” fungal phyla in both the P and N + P treatments when compared with the control and warming treatment groups (Figure 1b). In contrast, the relative abundance of the Basidiomycota phylum appeared greater in the birch rhizosphere of the N, P, and N + P treatments compared with the control and warming plots (Figure 1b, rhizosphere panel). Similar to the bacterial results, the average relative abundances of fungal phyla in the summer warming and control plots of the bulk soil appeared to be very similar. Supporting the observed differences in relative abundace of the fungi, ANCOM showed that Basidiomycota were differentially abundant in both the soil and rhizosphere treatments (Table S2). In addition, the Ascomycota and “other” phyla appeared differentially abundant among treatments in the rhizosphere (Table S2).

Changes in bulk soil and birch rhizosphere communities at the family and genus levels

Within both the bulk soil and rhizosphere samples, apparent differences in the average relative abundances of taxa at the family and genus levels were observed between treatments (Figure 2a). One of the more visible differences between treatments was a higher average relative abundance of Xanthomonadaceae (Figure 2a) in the N + P addition treatment plots, in particular when compared with the control and warming treatments in both the bulk soil and birch rhizosphere communities. Furthermore, this visible difference in relative abundance was noted in the separate N and P treatments compared with the N + P treatment. Fertilization (N, P, or N + P) treatments appeared to show a greater average relative abundance of soil and rhizosphere Sphingobacteriaceae compared with control and warming treatments and also resulted in a relative decrease in soil Planctomycetaceae when compared with the warming treatment (Figure 2a). N addition enhanced the average relative abundance of Chitinophagaceae in both soil and rhizosphere. In contrast, bacterial taxa from the orders GP1 and GP2 indicated a lower average relative abundance in all nutrient addition treatments (N, P, and N + P) compared with the control and warming plots for both soil and rhizosphere (Figure 2a). Similarly, at the phylum level, the bacterial bulk soil and rhizosphere community profiles in these plots appeared very similar in comparison to the other three fertilizer treatments (Figure 2a). Although there were visible differences in relative abundance between treatments, not all of these were supported by differential abundance analysis with ANCOM. However, several Acidobacteria-related families appeared to differ in relative abundance among the treatment groups (Table S3, Figure 2a). Additionally, the family with the highest differential abundance score is the Comammonadaceae in both the bulk soil and rhizosphere, which were present at <1 percent relative abundance across treatments (Table S3). The Comammonadaceae are higher in abundance in the N and N + P treatments, likely due to their involvement in N cycling.

Figure 2. Average relative abundances of taxa (family or genus) in bulk soil and birch rhizosphere within each experimental treatment group for both the (a) bacterial and (b) fungal communities of low Arctic mesic birch hummock tundra. Sequence results were analyzed to family level for (a) bacterial communities and genus level for (b) fungal communities in bulk soil and rhizosphere. The relative abundances shown represent the average of treatment replicates (n = 5), with treatment groups indicated on the x-axis. Any taxa present at less than 1 percent relative abundance across all treatment groups were pooled together for both the fungal and bacterial microbiomes.

Fertilization treatments also suggested differences in overall fungal community composition at the genus level compared to the associated controls, with N, P, and N + P treatment plots showing differences both in soil and the birch rhizosphere (Figure 2). For example, N + P amendment showed greater average relative abundance of the bulk soil genera Russula, Pseudogymnoascus, and Mycena and rhizosphere fungi Russula and Phialocephala (Figure 2b). P addition treatment indicated enhanced average relative abundance in the bulk soil of fungi belonging to the genera Thelephora, Mycena, Tomentella, and Ryvardenia and suggested a decrease in the relative abundance of Pezoloma (Figure 2b, soil panel). The P addition (and N + P treatment) also indicated an enhanced average relative abundance of the two fungal taxa Thelephora and Phialocephala in the rhizosphere (Figure 2b). With N amendment alone, Laccaria appeared to be relatively more apparent in the bulk soil, whereas both the N and P addition treatments suggested higher relative abundance of Mycena in the rhizosphere (Figure 2b). Once again, just as for the bacterial communities, the overall average relative abundances of fungal genus-level taxa between the warming and control treatments in either the bulk soil and or rhizosphere appeared similar, especially in contrast to the fertilizer-amended treatments. As with the bacterial family-level results, although some of the apparent differences in fungal genera between treatments were supported by ANCOM, others were not. The analysis did support that Thelephora, in both the bulk soil and rhizosphere, differed in abundance among the treatments and were particularly enhanced in the P treatment (Table S4). Likewise, Tomentella differed in abundance in the bulk soil, again due to substantial apparent increased abundance in the P treatment (Table S4).

Alpha diversity analysis

Faith’s PD as well as the Shannon index were used to determine the phylogenetic relationship between bacterial OTUs as well as the richness and evenness of bacterial taxa present within the samples. For bulk soil samples, the average alpha diversity in the bacterial microbiota indicated significant relative differences between treatment groups in both PD (p = .01) and Shannon diversity (p = .01), with an overall trend toward a decrease in diversity with P and N + P amendments (Figure 3a). The data suggested that among treatment groups, there was a small but significant relative decrease in PD in the N + P treatment plots compared with the controls (p = .025) and the N addition plots (p = .04). Similarly, the experiments indicated that there was a significant relative decrease in Shannon diversity in the N + P amendment compared with control (p = .046) and warming (p = .013) treatments. The data from the birch rhizosphere bacterial communities suggested that there were also significant relative differences in alpha diversity among treatment groups for PD (p = .02) and Shannon diversity (p = .01), with an apparent significant decrease in PD and Shannon diversity observed in the N + P treatment plots compared with the N treatments (p = .02) or controls (p = .02), respectively (Figure 3).

Figure 3. Phylogenetic and ecological diversity in bulk soil and birch rhizosphere microbial communities in arctic mesic birch soil for each experimental treatment group. (a) Average alpha diversity was calculated for bacterial diversity using PD and Shannon index—that is, the combination of richness plus evenness (SD)—and (b) fungal diversity was calculated using the Chao1 index and Shannon index (SD; n = 5). Asterisks indicate that the treatment group was significantly different from controls (p < .05). No significant relative abundance differences were found between treatments in fungal communities (p > .05) .

For fungal communities, alpha diversity was calculated using the Chao1 index but also with the Shannon diversity index, and no significant relative differences (p > .05) were seen between treatment groups in the bulk soil or the birch rhizosphere (Figure 3b). Although there was a significant p value noted for the Chao1 index when comparing soil fungal communities (p = .023), post hoc analysis indicated no significant differences between specific treatments.

Beta diversity analysis

The beta diversity analysis for bacterial microbiota and the Bray-Curtis metric for the fungal microbiota were plotted using NMDS to more readily visualize relationships between the treatment groups. The NMDS plots demonstrated a clear clustering of both the bulk soil and birch rhizosphere bacterial communities according to each treatment group (Figure 4b). The statistical significance of these cluster patterns was confirmed using nonparametric tests for similarities (ANOSIM) as well as for differences (adonis; Table 2; p < .001). It is notable that control and warming treatment groups overlapped extensively and were statistically grouped together (Tables S5 and S6; p > .05) and that all of the fertilizer treatment plots were grouped elsewhere, with the N + P amendment plots in particular located in a distinctly different cluster. Beta diversity distances within and between treatments indicated that soils from the warming and control treatments showed significant relative differences from both the N and N + P treatment groups (Tables S5 and S6; p < .05). Similar clustering patterns and corresponding statistically significant differences were also observed in the birch rhizosphere bacterial samples (Figure 4b).

Table 2. Analysis of similarities (ANOSIM) and of differences (adonis) between the bulk soil and birch rhizosphere beta diversity values for the bacterial and fungal communities across all experimental and control plots described in this study

Community	Type	Test	Test statistic*	p Value
Bacterial	Bulk soil	ANOSIM	0.651	.001
		adonis	0.640	.001
	Birch rhizosphere	ANOSIM	0.710	.001
		adonis	0.624	.001
Fungal	Bulk soil	ANOSIM	0.565	.001
		adonis	0.418	.001
	Birch rhizosphere	ANOSIM	0.589	.001
		adonis	0.402	.001

Note. *Test statistics are R for ANOSIM and R² for adonis tests.

Figure 4. Beta diversity of samples is shown with NMDS plotted with fitted environmental variables for bacterial taxa in (a) bulk soil and and (b) birch rhizosphere communities using the weighted Unifrac metric and fungal taxa in (c) bulk soil and (d) birch rhizosphere using the Bray-Curtis metric. Environmental vectors are pointed in the direction of greatest change, with the ellipses representing the standard error from the centroid of each treatment group cluster. Treatment groups are indicated by the legend at the bottom, with environmental variables including nutrient concentration and pH as shown. Note that all samples were CSS normalized prior to beta diversity calculations. Red arrows indicate significant correlation of environmental variables (p < .05) with beta diversity when fitted using vegan (R package) .

Mantel correlation analyses between beta diversity and soil concentrations of NO₃-N, NH₄-N, PO₄-P, and pH showed that soil nutrient availability appeared to drive clustering among treatment groups (Table 3). Correlations with other soil factors such as dissolved organic N (DON), microbial biomass N (MBN), dissolved total nitrogen (DTN), microbial biomass C (MBC), and microbial biomass P (MBP) were also performed. The data indicated that bacterial beta diversity in the bulk soils was significantly correlated with NO₃-N and PO₄-P concentrations as well as pH (p < .05). Furthermore, increases in NO₃-N and PO₄-P appeared to positively correspond with the direction of the clustered N and P treatment groups, respectively (Figure 4a). Significant correlations were also found with DON, DTN, and MBC (p < .05). In the birch rhizosphere samples, there was a significant correlation between bacterial beta diversity and PO₄-P concentration (p < .05), which correlated with the clustering of P treatment samples (Figure 4b). Although there appears to be some correlation between the clustering of the rhizosphere samples from the N treatment and soil pH, NO₃-N, and NH₄-N, as indicated by the vector arrows (Figure 4b), Mantel tests did not show any significant correlation between these factors and beta diversity across all treatments (Table 3), suggesting that there may be some selective process acting on the rhizosphere communities. The relative beta diversity of bacteria in the rhizosphere was also significantly correlated with MBN, MBC, and DTN (p < .05).

Table 3. Mantel tests results for only the significant correlations between beta diversity and biogeochemical properties at Daring Lake experimental plots in this study

Community	Sample type	Factor	R statistic	p Value^a
Bacterial	Bulk soil	MBC	0.272	.002
		DTN	0.214	.012
		NO3	0.248	.012
		pH	0.182	.034
		DON	0.174	.041
		PO4	0.176	.041
	Birch rhizosphere	MBC	0.325	.001
		PO4	0.244	.005
		DTN	0.173	.020
		MBN	0.173	.031
Fungal	Bulk soil	MBC	0.269	.004
		NO3	0.277	.008
		PO4	0.316	.009
		DTN	0.200	.031
	Birch rhizosphere	MBC	0.420	.001
		PO4	0.410	.003

Note. ^aFor both the bacterial and fungal soil and rhizosphere, only significantly correlated factors are included in the table (p < .05).

Beta diversity of the fungal community also suggested distinct clustering between some treatment groups (Figures 4c and 4d). This clustering was confirmed using adonis and ANOSIM analysis (p < .001). Again, the control and warming treatment groups overlapped with each other. This pattern was observed in both the bulk soil and rhizosphere communities where there appeared to be no significant difference in relative beta diversity between the two groups (p > .05). For the bulk soil, all three nutrient addition treatments (P, N, and N + P) visibly overlapped (Figure 4c), whereas the rhizosphere fungal microbiota showed that the N, P, and N + P treatments were clustered closely but not entirely overlapping (Figure 4d). Keeping in mind the experimental design, analysis of beta diversity distances within and between treatments confirmed that there was significantly more distance between the control samples and the N, P, and N + P treatments than within the control treatment (p < .05). The same was observed after comparison of the warming treatment to the three different nutrient addition treatments (Tables S5 and S6). There was, however, no significant statistical clustering among the nutrient-amended plots, as observed with beta diversity ordinations (Figures 4c and 4d). Again, Mantel testing using beta diversity and soil factors showed significant correlations with NO₃-N and PO₄-P, changing in the direction of the N and P treatments, respectively (Table 3). There was also correlation with MBC and DTN for the soil fungi, whereas the fungal rhizosphere community showed significant correlations only between MBC and PO_4, correlated in the direction of the P treatment cluster (Table 3).

Discussion

Long-term warming plots as simulations of climate change

Elevated temperatures associated with climate change are anticipated to gradually increase soil microbial activity in arctic terrestrial ecosystems and thus have significant impacts on vegetation. To better predict future changes, several long-term warming experiments have been established to simulate moderate climate shift. At the Daring Lake site, accumulation of plant aboveground nutrients increased in response to the summer greenhouse treatment (Zamin 2013). This suggests that warming likely increases decomposer microbial activity, resulting in the breakdown of soil organic matter. That soil microbial nutrient and biomass C pools were not altered by warming (Zamin 2013; Zamin, Bret-Harte, and Grogan 2014) implies that microbially mediated nutrient fluxes from soil organic matter into plant-available soil solution pools were enhanced without increasing microbial accumulation/retention of those nutrients. Here, our data suggest that the ongoing summer greenhouse warming regime did not significantly alter the relative compositions of the bulk soil or birch rhizosphere microbial communities compared to control plots, as indicated by overlapping clusters of taxa in both the soil and the rhizosphere of the warming and control treatment groups. Thus, observed increases in nutrient supplies to plants in warmed tundra soils (Zamin 2013) were not mediated by shifts in relative abundance of bacterial or fungal community members, with the most likely explanation being that warming enhanced biogeochemical activity of resident ambient microbial community members.

Why is it, then, that nutrients released as a result of such suggested temperature-mediated increases in soil microbial activity did not lead to obvious shifts in the soil community composition? Some long-term warming studies on microbial communities have shown alterations of both soil bacterial and fungal communities (Deslippe et al. 2012). However, other studies have demonstrated a decrease in soil microbial growth or no significant changes (Lamb et al. 2011; Rinnan, Michelsen, and Bååth 2011; Ballhausen et al. 2020), with still others showing increased microbial nutrient mobilization only if competing plants are eliminated (Schmidt et al. 2002). It is likely that both abiotic and biotic factors play a role in explaining these variable responses. For example, the Toolik Lake Alaskan warming treatment was conducted over two decades in moist soils and perhaps as a consequence showed more modest changes in mean summer experimental air temperatures compared to ambient controls relative to our experiment (mean increases of 1.3°C vs. 2.5°C; Deslippe et al. 2012 vs. Zamin 2013). In addition, compared to the Daring Lake site, the Toolik Lake site is more than 400 km further north and relatively wet with a deeper soil organic layer and characterized by vegetation that is a little different in that it is more dominated by tall Eriophorum tussocks. This may be of particular importance because previous results from our site indicated that the plant community (especially Rhododendron subarcticum and to a lesser extent B. glandulosa shrubs) outcompeted the soil microbial community for the modest levels of enhanced nutrient release in the greenhouse plots (Zamin 2013). Indeed, some arctic plants seem to take up nutrients more readily than others (Churchland et al. 2010), with certain species benefiting more than soil microbes from the enhanced N and P made available by warming (Schmidt et al. 2002). We found no evidence of significant warming-induced changes in the bulk soil or in the birch rhizosphere bacterial or fungal communities, suggesting that the compositional structure of the soil microbiota in our birch hummock tundra ecosystem was resilient to the greenhouse-mediated temperature increases. Overall, our data, in combination with previous studies (Zamin 2013), lead us to conclude that the greenhouse treatment increased soil organic matter microbial decomposition activity, thereby enhancing soil nutrient availability, which was subsequently accumulated primarily by plants rather than microbes. As a result, there was no significant alteration in the compositional structure of either the bulk soil or birch rhizosphere bacterial or fungal communities.

It would be interesting to continue monitoring these warmed plots in the future because significant impact of warming on soil communities may require longer experimental treatment time than a dozen years, although it is noted that not all researchers support longer-term experiments (Ballhausen et al. 2020). Alternatively, because temperatures in the greenhouses on rare extremely warm days can rise as much as 13°C above ambient (Zamin and Grogan 2012), combined with rising background climate change–associated increases, future warming experiments could provide unrealistic conditions for arctic soil communities and unduly stress the plants. This then suggests the need for some protective measures within the greenhouses to prevent such overheating or, alternatively, the use of open-top chambers.

Long-term nutrient additions as simulations of climate change

As is well understood, climate warming results in increased soil nutrient availability. Chronic annual large fertilizer addition experiments may simulate extreme warming effects but are generally interpreted as being most useful to identify the impacts of completely removing nutrient limitation on plant and microbial communities. At the Toolik Alaska site, such regimes have resulted in shifts from graminoid to shrub vegetation (Nowinski et al. 2008; Prager et al. 2017) and changes in soil bacterial structure but not in the fungal communities (Campbell et al. 2010; Koyama et al. 2014). In contrast, at the Daring Lake site both bacterial and fungal communities appeared to be substantially altered by direct nutrient amendments. The most striking changes were often seen with N + P additions showing a relative decrease in richness and evenness (alpha diversity; SD) in bacterial communities, as well as relative shifts in fungal abundances in the P and N + P amended treatment plots. It could be argued that different DNA primers could explain some of the variation in responses between the Canadian and Alaskan sites (Liu et al. 2015), but it is more likely that environmental and vegetation differences at the two sites, as previously discussed, play an important role, including water drainage, which is known to have a pronounced influence on C degradation and accumulation (Schuur et al. 2009; Dahl et al. 2017). We thus posit that the differences in fungal community responses between the two study locations can be attributed to site differences including vegetation as well as water potential due to flooding frequency and soil moisture levels, which, as noted, are higher at the Alaskan site through much of the growing season.

As previously mentioned and supported by a rich literature, abiotic factors are correlated with changes in bacterial and fungal structure (e.g., Malard and Pearce 2018). The beta diversity of bacterial soil communities is significantly impacted by soil pH, but here the lack of a correlation for this abiotic factor and the relative abundance of the microbiota of the rhizosphere communities likely reflects the increased consistent pH of the rhizosphere compared to bulk soil (Bagayoko et al. 2000; J.-Y. Shi et al. 2011). Nutrient additions resulted in relative increases in C, N, and P pools, which in turn influenced overall community profies (beta diversity), and both bacterial and fungal communities in the bulk soil and birch rhizosphere significantly correlated with MBC. The significant decrease in microbial biomass seen in the P and N + P amended plots (Zamin, Bret-Harte, and Grogan 2014) and the observed relative shifts in community structure suggest that the nutrient amendments reduced both community diversity and microbial biomass. Microbiota change has been linked to N availability (e.g., Ramirez, Craine, and Fierer 2012; Koyama et al. 2014; N. F. Wang et al. 2016). Thus, it is notable that in our plots N amendment did not appear to have as striking an impact on microbial structure as other nutrient additions. For example, NH₄-N concentrations were not significantly correlated with beta diversity, although NO₃-N was significantly correlated with soil bacterial and fungal communities. P concentrations were significantly correlated with beta diversity in each of the bulk soil and birch rhizosphere clusters for both bacterial and fungal microbiota. Thus, the availability of P, either by itself or with additional N, is an important factor driving the structure of both the bacterial and fungal communities in the bulk soil and in the birch rhizosphere of this arctic ecosystem. Indeed, the availability of P, either by itself or with N, is also known to be an important determinant of birch shrub growth at the Daring Lake site (Zamin and Grogan 2012; Zamin, Bret-Harte, and Grogan 2014).

Bacterial communities in both the bulk soil and birch rhizosphere samples from across our experiment appeared to be dominated by Actinobacteria and Proteobacteria, respectively, and similar patterns have been reported from other low Arctic as well as high Arctic and alpine tundra sites (Y. Shi et al. 2015; Kumar et al. 2017; Voříšková, Elberling, and Priemé 2019). Fertilization, particularly with N + P, seemed to decrease the average relative abundance of Planctomycetes, a taxon associated with lichens (Ivanova et al. 2016), and thus is consistent with the potential for landscape change. However, it should be noted that this observation was not confirmed using differential abundance analyses (Table S1), which are not suitable for post hoc analysis such as the comparison of pairs of treatment groups. Soil amendments with P and N + P increased the average relative abundance of Proteobacteria, such as Sphingobacteriaceae and Chitinophagaceae, suggesting again that P is likely important and may be a limiting nutrient for members of these two families. The apparent increase in the relative abundance of Xanthomonadaceae in the bulk soil and birch rhizosphere of N + P amended plots, predominantly of the genus Rhodanobacter, is consistent with the prevaence of these bacteria in other long-term fertilized soils and the association with arcto-alpine plants (Jangid et al. 2008; Campbell et al. 2010; Nissinen, Männistö, and van Elsas 2012; Li et al. 2017). Differential abundance analyses and these observations together further highlight the potential importance of these bacteria, including some denitrifying members within communities of nutrient-manipulated soils, and their possible utility as an indicator for relatively fertile arctic soils.

Our experiments suggest that fungal community structure could also be altered by fertilizer amendments. For example, within the P treatment groups there appeared to be a significant differential abundance in Basidiomycota, supported by ANCOM, with the biggest change being a relative increase in the ectomycorrhizal fungus Thelephora in both the bulk soil and birch rhizosphere samples. This result is consistent with the increase in this taxon after P amendments on pine trees, Pinus sylvestris (Colpaert et al. 1997). The data further suggest that within the N + P amendment plots, the largest average relative increase in fungal taxa compared to the associated controls was for members of Russula, another ectomycorrhizal fungus, which is a common associate of Arctic birch roots. Although not supported by ANCOM, there also appeared to be a relative increase in the Ascomycota Phialocephala after P and N + P amendments in the rhizosphere. Because Phialocephala iis a potential decomposer of soil organic compounds (Surono and Narisawa 2017), an increase in this endophyte could result in an increase in available C, which in turn would promote microbial activity, including ectomycorrhizal fungi and rhizosphere-associated bacteria, and, by extension, promote plant growth (Newsham 2011). We thus posit that fertilizer-mediated changes in microbiota could foster relative increases in fungal biomass or, alternatively, could stimulate the growth of ectomycorrhizal and rhizosphere fungi associated with B. glandulosa, which in turn could increase the fitness of this plant species, making it more competitive. Previously, mycorrhizal networks, which promote the sharing of resources among plants and fungi, have been associated with an increased competitive advantage for a number of plants (Simard et al. 2015). Furthermore, B. glandulosa is the only ectomycorrhizal host plant species in most low Arctic tundra ecosystems and certainly in the mesic birch hummock tundra at our site. Indeed, a strong increase in aboveground biomass of Arctic birch in these fertilization treatment plots has been reported (Zamin, Bret-Harte, and Grogan 2014), and here we propose that this likely reflects the observed changes in the fungal community. Arctic birch must be highly competitive and more growth responsive compared to many other tundra plant species, because it is currently increasing in density and expanding across regions of the Arctic (Chapin et al. 1995; Tremblay, Lévesque, and Boudreau 2012). Therefore, our fertilization treatment results suggest that in locations where tundra soil nutrient availability is markedly enhanced by either nutrient inputs such as caribou mortality or pollution or even possibly in some regions by warming, such changes in the fungal consortia could contribute to Arctic birch expansion.

Conclusions

Though other studies have evaluated the impact of long-term fertilization on arctic soils using combined N + P nutrient amendments (e.g., Campbell et al. 2010; Koyama et al. 2014), to our knowledge our experiments are the first to evaluate the individual impacts of N vs. P and compare them to N + P additions on microbial communities in the tundra. As previously noted, it is commonly regarded that N availability is an important driver of microbial diversity. However, we observed greater relative shifts with P and N + P amendments, compared with N only, suggesting that P may be a critical, although less acknowledged, driver of change in microbial communities. Certainly, our experiments indicate that nutrient addition significantly changes both bulk soil and rhizosphere communities, the latter of which can have important implications for health and growth promotion for Arctic birch, a vital component of arctic tundra ecosystems.

Overall our results also demonstrate the complexity of using experimental manipulations to try to understand how climate change will affect bacterial and fungal communities in tundra soils across arctic ecosystems. The absence of statistically significant greenhouse warming effects on soil and rhizosphere community composition, compared to the obvious relative differences observed in the nutrient-amended plots, clearly shows that these two types of treatments have profoundly different impacts. The former is constrained by factors such as artificially enhanced temperature amplitudes (especially maximums) in a single season of each year, even though the mean temperature increase reasonably simulated anticipated climate warming. The fertilization treatments increase nutrient availability to amounts that are higher than levels predicted from climate change over this century and instead may help us to understand and compare nutrient limitations on microbiota and plant growth. Both treatment types therefore allow an appreciation of the mechanistic factors determining the structure (i.e., richness, relative abundance, composition, and interactions) of tundra microbial and plant communities. Therefore, despite the caveats associated with such experimental treatments, we suggest that these manipulations, especially the low-level nutrient addition plots, will continue to provide insights into the mechanisms determining tundra microbial–vegetation communities and help us understand the consequences of climate warming on these ecosystems.

Acknowledgments

We acknowlege the excellence of the Centre for the Analysis of Genome Evolution & Function (CAGEF) at the University of Toronto and, in particular, Dr. D. Guttman and Julia Copeland. We thank Qian Gu for assistance with the fieldwork and K. Moniz for her valuable technical expertise. We acknowledge anonymous reviewers for their comments.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Funding

This research was supported by funding from the Natural Sciences and Engineering Research Council (NSERC) of Canada Discovery grants to V.K.W. and P.G.

Supplementary material

Supplemental material for this article can be accessed on the publisher’s website.