Consistent responses of soil microbial communities to elevated nutrient inputs in grasslands across the globe

Fungal diversity and community composition differed strongly across the 25 globally distributed grassland sites, regardless of nutrient treatment (P < 0.001 in all cases; Fig. S1). Mean fungal phylotype (i.e., species) richness ranged 1.7-fold across the sites, and there were large variations in the relative abundances of major taxonomic groups (Table S1). The strong site effects are not surprising, given the range in environmental conditions and soil characteristics found across sites spanning four continents and elevations from 50 to 2,320 m (Tables S2 and S3). In particular, the sites represented a broad range in soil acidity, climate, and plant community composition, factors that have previously been associated with differences in soil fungal community structure at these sites and others (30, 31).

We investigated the within-site effects of nutrient additions on fungal community structure by statistically controlling for the strong cross-site differences by including site as a random effect in our models. Fungal Shannon diversity responded weakly to nutrient additions, decreasing by only 2.7% on average when N and P were added together (P = 0.05), a response consistent with the weak response observed for plants (8).

In contrast to the weak effects of nutrients on fungal diversity, we observed significant effects of N (R² = 0.003; P < 0.001) and/or P (R² = 0.002; P = 0.04) additions on fungal community composition, with the same taxa generally responding to nutrient additions across sites, despite the large cross-site variation in fungal community types (Fig. 1). With combined addition of N and P, there were increases in Ascomycota and significant decreases in the relative abundances of Glomeromycota (Fig. 2A). The Glomeromycota phylum is composed almost entirely of arbuscular mycorrhizal fungi (32), and we expected these fungi to decrease in relative abundance with nutrient additions because they would be less valuable to their hosts and thus provided with less plant C under conditions of increased N and P availability (33–35). We further investigated nutrient effects on mycorrhizal fungi by assessing the collective responses of mycorrhizal fungi, including those taxa outside the Glomeromycota phylum that are reported in the literature as being mycorrhizal. These taxa also consistently decreased in plots receiving N and P relative to the control plots (P = 0.016), corroborating results from a meta-analysis demonstrating declines in mycorrhizal fungi with N additions (18). Interestingly, adding N and P together led to far larger decreases in the relative abundances of Glomeromycota than when these nutrients were added individually (P > 0.1; Table S4), suggesting a role for both of these nutrients in shaping arbuscular mycorrhizal communities.

The overall decrease in the proportion of mycorrhizal fungi with N and P additions—and shifts in fungal community composition more broadly—could be caused by plant community shifts, changes in plant biomass, and/or the direct effects of added nutrients. The magnitudes of the responses of major fungal taxonomic groups were not significantly correlated with changes in key soil characteristics (Table S5). However, the magnitude of fungal community composition response (i.e., the mean community dissimilarity between samples with added N and P and control samples) was significantly correlated with the magnitude of the response of plant community composition to added N and P (r = 0.44; P = 0.03; Fig. 3), helping to explain site-to-site variability in shifts in belowground communities. Those sites where nutrients had the largest impacts on plant communities were also the sites that had the strongest nutrient effects on fungal communities. This finding suggests either that shifts in plant community composition drive shifts in fungal community composition or that both plant and fungal communities respond similarly to changes in edaphic factors. Although overall fungal compositional shifts correlated with plant community composition shifts, changes in the relative abundance of Glomeromycota were not related to changes in live plant biomass with fertilization (P > 0.1), or to changes in surface soil N concentrations (P > 0.1; Table S5), suggesting that plant nutrient limitation was not a good predictor of the differential responses observed across the sites.

Archaea were rare at most sites, and archaeal diversity (Fig. S1A) and community composition (Fig. S1B) were highly variable across sites, regardless of nutrient additions (P < 0.001). Archaeal phylotype richness ranged 3.7-fold across the sites, and the archaeal communities were dominated by Crenarchaeota (92% on average) and Euryarchaeota (4.3% on average; Table S1). The proportion of 16S rRNA reads that were of archaeal origin was also highly variable across the sites (Fig. S2A), ranging from 0 to 0.16. This variability in archaeal communities was likely due to the large cross-site differences in environmental conditions mentioned above. For instance, previous work has shown a correlation between archaeal relative abundances and soil nutrient content (36); we know that soil N concentrations varied 33-fold across the control plots, and archaea relative abundances were inversely related to soil C:N ratios (r = −0.67; P < 0.001).

We next assessed whether there were consistent shifts in archaeal relative abundance and community structure with nutrient additions by statistically controlling for the strong cross-site differences. Archaeal relative abundances generally increased with N additions (P < 0.001; Fig. S2B), and there was a mean 4.8% decrease in archaeal diversity with N additions compared with control plots (P = 0.01). This decrease in diversity was possibly related to an N-induced growth of specific archaeal taxa. Specifically, the phylum Crenarchaeota, which primarily comprised members of the family Nitrososphaeraceae, consistently increased in relative abundance with N additions across the majority of sites, whereas Euryarchaeota and the candidate division Parvarchaeota consistently decreased (Fig. 2B). These shifts are likely related to Archaea being active drivers of the soil N cycle. For example, Nitrososphaeraceae can oxidize ammonia (29, 37), a metabolism that is expected to be advantageous with elevated ammonium supply, which should have been elevated in the N addition plots, because urea is readily hydrolyzed to ammonium. Abundances of soil Crenarchaeota also are positively correlated with soil N content (36). Conversely, several reports have shown the potential for members of the Euryarchaeota, which are predominately methanogens, to fix atmospheric N₂ (38, 39). This characteristic could place them at a competitive disadvantage under conditions of elevated N availability and explain their strong proportional decrease with N fertilization. Although it has been shown that N can inhibit methanogenesis in vitro (40), this work is, to our knowledge, the first direct evidence that N additions may also decrease methanogen populations in non-wetland soils. Still, it is important to note that these shifts in the relative abundances of archaeal phyla are not independent of one another, and decreased methanogen relative abundances could simply be the result of increased relative abundances of Crenarchaeota. Nonetheless, these results highlight that soil archaeal communities are sensitive to N additions, but additional research is required to determine whether these community responses are associated with changes in methane fluxes or soil N cycling rates.

As with fungal and archaeal communities, bacterial diversity and community composition differed strongly across the 25 grassland sites (Fig. S1). These differences were likely due to factors such as acidity, climate, and plant community composition, as has been previously observed (30, 41, 42). Mean phylotype richness ranged 1.7-fold, and the abundant phyla, including Proteobacteria, Acidobacteria, Verrucomicrobia, Actinobacteria, and Bacteroidetes, all varied considerably in their relative abundances across the sites (Table S1).

Nutrient additions did not strongly alter bacterial diversity; P additions caused marginal (0.5%) increases in bacterial diversity (P = 0.06), and N had no significant effect. Our results stand in contrast to negative relationships between bacterial diversity and N additions reported from previous studies conducted at individual sites (19, 43). This finding points to the importance of local context and highlights the pitfalls associated with extrapolating results obtained from individual sites to other ecosystems or soil types.

Bacterial community composition was significantly affected by N (R² = 0.002; P < 0.001) and P additions (R² = 0.002; P = 0.003; Fig. 1). The community shifts corresponded to changes in the relative abundances of numerous major taxa. For example, the relative abundances of Actinobacteria, Alphaproteobacteria, and Gammaproteobacteria consistently increased with nutrient additions across sites, whereas those of Acidobacteria, Planctomycetes, and Deltaproteobacteria consistently decreased (Fig. 2C). However, these taxonomic shifts were not always in the same direction or magnitude when N or P was added alone (Table S4). Overall, the taxonomic patterns in our cross-site study were in agreement with previous work conducted at individual grassland sites (20), and they corroborate laboratory studies that have noted similar shifts in the relative abundances of these major bacterial groups with nutrient additions (13). Our findings are generally consistent with our hypothesized shifts in general life history strategies with bacterial taxa that are faster growing and more copiotrophic (28) being favored under conditions of elevated nutrient availability (27). In particular, soil bacterial groups that are generally considered to be more copiotrophic, including Actinobacteria and Alphaproteobacteria, increased in relative abundance with nutrient additions, and the largely oligotrophic Acidobacteria phylum decreased in relative abundance. Whereas original evidence for generalizations of these life history strategies across broad bacterial taxonomic groups was based on responses to labile C inputs (28, 44, 45), our results extend evidence for these ecological classifications to the direct or indirect bacterial responses to nutrient additions.

We recognize that it is difficult to confidently assign bacterial clades into groups with copiotrophic and oligotrophic life history strategies, especially given the overwhelming amount of undescribed bacterial diversity found in soil (25). Thus, we used a combination of genomic and metagenomic approaches to provide independent assessments of how copiotroph:oligotroph ratios shifted in response to added nutrients. First, we estimated aggregate community growth rates because we expected increases in the relative abundance of copiotrophic taxa to be reflected by faster growth rates (28, 46). Thus, an increase in the estimated growth rate [i.e., a decrease in mean minimum generation time (MGT)] would suggest an increase in the relative abundance of copiotrophs. Mean MGTs were calculated for all samples from a combination of our bacterial marker gene data and published genomes; 757 of the 46,534 phylotypes could be matched to genomes. As with other attributes of community structure, estimates of MGT strongly varied across sites (Fig. S3A). Within-site differences between nutrient-amended and control samples showed that adding nutrients tended to decrease MGTs (Fig. S3B), but this trend was not significant for N additions (P = 0.57) or P additions (P = 0.34) individually. However, this analysis has important limitations in that only a small proportion (∼10%) of the 16S rRNA gene sequences from our samples could be mapped to genomes for which we had MGT estimates, and this proportion differed across nutrient treatments (Fig. S3C). Thus, this analysis likely provides a conservative estimate of potential differences in MGTs associated with nutrient additions and is only weakly supportive of the hypothesis that soil bacterial MGT decreases with nutrient additions.

To further confirm the putative shifts in life history strategies in bacterial communities, we assessed functional attributes directly from functional gene (i.e., shotgun metagenomic) data collected from six of the sites used in the taxonomic analyses (Tables S2 and S3). These sites were selected because they spanned a wide geographic range and encapsulated a variety of environmental conditions, and the marker gene analyses suggested the N and P effects on microbial community composition were particularly strong. The shotgun metagenomic data (hereafter referred to as metagenomic data) were found to be almost entirely derived from bacterial genomes—94.8 ± 2.3% (mean ± SD) of the metagenomic small subunit rRNA gene reads were identified as bacterial. Just as the marker gene data revealed that bacterial diversity and community composition differed strongly across sites, the metagenomic data revealed that functional gene diversity and composition also varied strongly across sites (Fig. S1). In addition, the diversity of annotated genes identified from the metagenomic data were significantly correlated with the diversity of bacterial phylotypes across the samples (r² = 0.27, P < 0.001; Fig. S4A), and the dissimilarity in functional gene composition was strongly related to the dissimilarity in bacterial community composition across the six sites (ρ = 0.87, P < 0.001; Fig. S4B). These findings suggest that bacterial communities that are distinct in composition tend to have distinct functional attributes, and bacterial communities that are taxonomically more diverse also have more diverse metagenomes with a broader array of annotated genes. Correspondingly, the diversity of functional genes did not change with nutrient additions (P > 0.1), but there were significant shifts in overall functional gene composition with N additions (P = 0.01) and P additions (P = 0.006; Fig. 1), as was observed for bacterial taxa. These results are supported by previous work showing a relationship between the taxonomic structure of soil bacteria and functional genes across ecosystems (41) and significant N effects on functional gene composition at two North American sites (27).

The metagenomic data yielded additional lines of evidence to support our hypothesis that nutrient additions favor copiotrophic bacterial taxa. Previous work has suggested that soil microorganisms with larger genomes should be more successful in resource-poor environments (47), and, thus, we expect copiotrophic taxa to have smaller genomes. To assess this hypothesis, we calculated mean effective genome size—the estimated mean size of a genome in a given sample—and found that it significantly decreased with added N or P (P < 0.03 in both cases; Fig. 4A). More generally, this result highlights that genome size can be considered an important ecological trait, just as bacterial genome size is correlated with range size (48) and plant genome size is an important predictor of species’ ability to invade (49).

We investigated the specific gene categories that changed in proportion with nutrient additions by analyzing the quality-filtered metagenomic sequences that could be annotated. First, it is important to note that only 28.7–32.7% of sequences could be annotated, and soils receiving N or P had a 0.3% higher annotation rate on average (P ≤ 0.01 in both cases; Fig. 4B), a pattern likely driven by the overrepresentation of copiotrophic bacteria, which are easier to culture, and are thus more commonly found in genome databases. Similarly, soils receiving N amendments tended to have a lower relative abundance of annotated, but unclassified, metabolic genes compared with control samples, likely also reflecting the better representation of copiotrophs in genome databases (Fig. 4C and Table S6). We also observed a significant increase in the relative abundances of genes associated with carbohydrate metabolism (Fig. 4C) in fertilized plots. This finding is consistent with the added nutrients increasing copiotroph:oligotroph ratios and potentially increasing plant C inputs to soil. Although <33% of the sequence reads could be annotated—a percentage that is similar to that reported in other metagenomic analyses of diverse bacterial communities, e.g., ref. 27—our results highlight that the annotated reads can be used to infer shifts in the functional capabilities of communities, shifts that are consistent with nutrient additions increasing the proportional abundance of bacteria with copiotrophic life history strategies.

Nutrients can have both direct and indirect effects on belowground bacterial communities, making it difficult to unravel the mechanisms underlying the community responses described above. Potential mechanisms include direct effects of the nutrients themselves, nutrient effects on soil characteristics (e.g., pH), nutrient inputs increasing plant productivity and organic matter inputs to soils (20), and nutrient inputs mediating microbial shifts through changes in plant community composition. With N addition, soil pH decreased by an average of 0.16 units across the sites (P < 0.001), and pH has been shown to strongly drive shifts in soil bacterial communities (42, 50, 51). However, pH alone is not likely to have been a major driver of community shifts observed here, because the pH change was relatively small, it did not change with P additions (P = 0.36), and the magnitude of change in pH was unrelated to the change in the relative abundance of any of the major bacterial taxa with N and P additions across the sites (Table S7). Proportional changes in plant productivity were also unrelated to changes in the relative abundance of bacterial taxa, suggesting that elevated plant productivity in fertilized plots was not responsible for the bacterial community responses. Conversely, the magnitude of shifts in plant community composition was directly related to the magnitude of shifts in bacterial community composition (r = 0.41, P = 0.04; Fig. 3), a pattern that mirrored that observed for fungi (Fig. 3). These findings suggest that changes in plant community composition may be more important for mediating bacterial community responses to elevated nutrient inputs than changes in edaphic characteristics or plant growth.

All 25 sites used in this study (Tables S2 and S3) were from temperate-zone grasslands in Africa, Australia, Europe, and North America, and all were part of the Nutrient Network experiment (8) and are described in ref. 30. In brief, the sites ranged in many environmental characteristics, including mean annual precipitation (262–1,898 mm⋅y⁻¹), mean annual temperature (0–18 °C), elevation (50–2,320 m), soil pH (4.5–8.4), total soil P (1–253 ppm), soil %N (0.03–1.5%), and aboveground plant productivity (15–1482 g⁻²⋅y⁻¹).

Identical full factorial N and P addition experiments were established at each of the 25 sites (8). We used samples from three +N, three +P, and three +N/+P 1-m² plots, and six 1-m² plots with no added nutrients. Nutrients were added annually as 10 g of N per m²⋅yr⁻¹ timed-release urea [(NH₂)₂CO] and 10 g of P per m²⋅yr⁻¹ triple-super phosphate [Ca(H₂PO₄)₂]. Samples were treated between 2 and 4 y before collection (Table S2), but this timing did not significantly relate to the proportional changes in the major microbial taxon relative abundances (P > 0.5 in all cases). Fences to exclude herbivores surrounded three of the plots at each site (8), but because fencing had no effect on belowground microbial communities (P > 0.1), we included these plots as controls.

Plant community composition and biomass were assessed and soil samples were collected during the growing season of 2011 or 2012. Plant species were identified within the plots, and plant biomass and soil cores were sampled directly adjacent to the plots to determine soil characteristics and assess microbial community structure (8). Soil cores were shipped on ice to a central processing facility (Corvallis, OR) immediately after collection, and samples for microbial community analysis were preserved at −20 °C. Analyses to determine soil pH and C, N, and P content are described in ref. 30.

Microbial community diversity and composition were assessed by using targeted marker gene surveys focusing on the 16S rRNA gene for the Bacteria and Archaea domains and the internal transcribed spacer (ITS) region for Fungi as in Prober et al. (30). DNA was extracted by inserting a sterile swab into each soil sample and cutting the swab tip off into a well in a bead plate of the PowerSoil-htp 96-well DNA extraction kit (Mo Bio Laboratories, Inc.). The downstream extraction procedure was followed according to the manufacturer’s instructions, and two wells on each extraction plate were left empty to serve as blanks. Marker genes in isolated DNA were PCR-amplified and barcoded in triplicate reactions for both the 16S rRNA gene (using the 515f/806r primer pair) and the ITS1 region (using the ITS1-F/ITS2 primer pair). PCR products from the triplicate reactions were combined and visualized on an agarose gel to ensure successful amplification and to verify no amplification from blanks. PCR product from samples that were successfully amplified was combined in equimolar ratios for each of the marker genes. This procedure was conducted separately for two sets of samples, one including samples from 19 of the sites and the other including the samples from the remaining 6 sites. For the first set, 16S rRNA and ITS amplicon pools were each sequenced on separate Illumina HiSeq 2000 lanes by using 100-bp paired-end sequencing, and for the second set, each amplicon type was sequenced on a separate Illumina MiSeq run by using 151-bp paired-end sequencing. All sequencing was conducted at the University of Colorado at Boulder.

Raw sequences from 16S rRNA and ITS amplicons were processed by using the UPARSE pipeline (58). Sequences were demultiplexed according to the raw sequenced barcodes, and because there were generally many more sequences for samples sequenced on the HiSeq run, a random subset of the sequences from each of these samples was included in further processing such that the number of sequences per sample was generally within two orders of magnitude across samples. Sequences from both sequencing runs were combined after trimming to a uniform length of 100 bp. A de novo database of ≥97% similar sequence clusters was created in USEARCH (Version 7; ref. 59) by (i) quality filtering sequences using a “maxee” value of 0.5 (i.e., sequences with a predicted error rate of 0.5 bases per sequence were discarded), (ii) dereplicating identical sequences, (iii) removing singleton sequences, (iv) clustering those sequences, and (v) filtering poor-quality sequences by removing sequences that were not ≥75% similar to any sequence in Greengenes (Version 13_5; ref. 60) or UNITE (Version 12_11; ref. 61) databases for 16S rRNA and ITS1 sequences, respectively. Raw demultiplexed sequences were then mapped against these de novo databases to generate counts of sequences matching clusters (i.e., phylotypes) for each sample. Taxonomy was assigned to each phylotype by using the RDP classifier with a confidence threshold of 0.5 (62) and trained on the databases indicated above. Among 16S rRNA sequences, phylotypes classified as bacteria were separated from those classified as archaea before further processing. Mycorrhizal fungi phylotypes were identified by comparing to known mycorrhizal ITS sequences, as detailed in ref. 30. To normalize the sequencing depth across samples, samples were rarefied to 18,000 bacterial, 100 archaeal, and 485 fungal sequences per sample. Because of insufficient sequence coverage, bacterial data from 4 samples, archaeal data from 90 samples, and fungal data from 29 samples were discarded. One sample (“NN10”) was removed from downstream analyses due to a large disparity in diversity and community composition from all other sequences. The raw sequence data are available in the Sequence Read Archive at the National Center for Biotechnology Information (BioProject accession no. PRJNA272747).

All replicate samples from a subset of six of the sites (90 samples; Fig. S1) used for the compositional analyses were prepared and submitted for shotgun metagenomic sequencing at the Joint Genome Institute. The sites were selected because they represented a broad range in environmental characteristics and exhibited the largest bacterial community composition responses to N and P. For this analysis, DNA was extracted by using the same method as above, except that ∼0.25 g of each sample was loaded into the first extraction plate. Library preparation and sequencing were conducted at the Department of Energy Joint Genome Institute (Walnut Creek, CA) using their standard protocols. Briefly, 270-bp DNA fragments were prepared and sequenced on a HiSeq instrument using 150-bp paired-end sequencing. Sequencing efforts yielded a total of 850 million read pairs across all samples (range = 428,000 to 65 million). Raw sequence data can be accessed at IMG (img.jgi.doe.gov) and are referenced in the Genomes Online Database (GOLD Study ID Gs0053063).

Sequence processing was conducted in a method similar to Fierer et al. (55). Paired-end sequences were first merged by using FLASH (Version 1.2.6; ref. 63) with the parameters, “-r 150 –f 210 –s 21.” Merged sequences were quality-filtered by using PRINSEQ-lite (Version 0.20.3; ref. 64) to remove sequences containing more than five undefined bases and exact duplicates. On average, 49% (range = 23–62%) of sequences were successfully merged and passed quality filtering from those samples that were successfully sequenced. Sequences from samples with high sequence coverage were then randomly subsampled to reduce computational processing time. The merged and quality-filtered sequences were annotated by mapping them against the IMG (Version 350) amino acid sequences (65) using BLAT (66). Sequences that could be mapped within a minimum percentage identity of 0.55 and an e-value cutoff of 1e-03 were provided Kyoto Encyclopedia of Genes and Genomes Orthology (KO) identifiers. Of the quality-filtered sequences, 28.7–32.7% could be annotated with specific gene category information, a similar percentage to that noted in previous studies exploring soil metagenomes (27), but larger than studies using comparatively short sequences (41). Only sequences that were annotated were considered for further analysis, and samples were rarefied to an equal sequence depth of 500,000 annotated sequences per sample. This process resulted in a total of 75 shotgun metagenomic samples spread across treatments and sites. When KO identifiers were categorized into coarser functional categories (56), those that belonged to multiple categories were counted in each of those categories. Functional diversity and composition calculations were performed identically to those used for taxonomic analyses, except that KO identifiers were used instead of phylotypes.

To determine what proportion of the shotgun data represented bacterial genes, we compared merged and quality-filtered reads to the SILVA 111 database (67), which contains representative small subunit rRNA sequences from all three domains of life. Shotgun reads were mapped to this database at the 94% similarity level by using UCLUST (59). We used a low similarity threshold to account for the diverse soil taxa that are not represented in the database.

To estimate mean MGTs (i.e., community-aggregated growth rates) from 16S rRNA sequence data, MGTs were estimated from published genomes and linked to bacterial community composition by matching 16S rRNA gene sequences. Fasta files of all ORFs from 1,136 genomes were obtained from the Joint Genome Institute’s IMG database. These genomes were selected to only include non-host-associated organisms based upon information available in the GOLD. Both IMG and GOLD were accessed in September 2012. MGT was estimated for each genome based upon the methods described by ref. 57 and python code provided by the authors. Briefly, two measures of codon use bias (ENC′ and S) were combined via regression to generate a predictive model of MGT. MGTs were linked to phylotypes observed through 16S rRNA data by matching the 16S rRNA sequences in our dataset to the 16S rRNA sequences from the genomes at a 97% similarity threshold using the “usearch_global” function in USEARCH (59), with “maxaccepts” and “maxrejects” set to 0 to disable approximate clustering. Mean MGT was calculated on a per-sample basis by calculating proportional abundances of the phylotypes that matched represented genomes, excluding all other phylotypes, and multiplying its estimated MGT by its proportional abundance. The mean MGT was the sum of these values.

Microbial diversity was calculated by using Shannon diversity because this method has been found to be a more reliable metric than species richness with microbial sequence data from complex communities (68). However, richness values are provided for interpretation and comparative purposes. Kruskal–Wallis tests were used to test whether there were significant differences in diversity among sites. Linear mixed effect models were used to test for differences in diversity with nutrient additions (+N, +P, +NP) by using Genstat (Version 14.0). We excluded the +NP term if not significant and treated site as a random variable.

We tested for significant shifts in overall community composition or functional gene composition across sites or with nutrient additions using permutational multivariate ANOVA (PERMANOVA) implemented in the “Adonis” function in the vegan package in R (Version 3.0.2; ref. 69). Community or genetic composition was represented by Bray–Curtis dissimilarity matrices computed from square-root-transformed abundance tables. N and P addition were used as predictor variables, and the geographic site was included as “strata,” which restricts permutations to within sites. N × P interactions were found to be nonsignificant in all cases and therefore removed from final PERMANOVA models. When visualizing potential treatment differences, we used constrained ordination as implemented in the “capscale” function in the vegan package in R. To test for significant differences in the relative abundance of specific taxa, gene categories, and MGTs between control and nutrient treatment plots, we used linear mixed-effects models in which each nutrient treatment (N and P additions) was a fixed factor and site was a random effect. Relative abundances were rank-transformed to meet the assumptions of the model. Tests were run for each taxon represented by a median of at least 1% of the sequences in any of the nutrient treatments, and false discovery rate (FDR) corrections were used to account for the multiple comparisons.

Relationships between pretreatment environmental variables and changes in the relative abundance of taxa were assessed by using Pearson correlations with mean difference from control plots (on a per-site basis) as the response variable. We chose environmental variables to test based on their relationships with one another and avoided including multiple variables that were highly correlated with one another (r² > 0.8). Differences in posttreatment soil pH were investigated by using a linear mixed-effects model with N and P additions as fixed effects and site as a random effect. This analysis was conducted by using the “nlme” package in R.