Chronic N-amended soils exhibit an altered bacterial community structure in Harvard Forest, MA, USA

Abstract

At the Harvard Forest, Petersham, MA, the impact of 20 years of annual ammonium nitrate application to the mixed hardwood stand on soil bacterial communities was studied using 16S rRNA genes pyrosequencing. Amplification of 16S rRNA genes was done using DNA extracted from 30 soil samples (three treatments × two horizons × five subplots) collected from untreated (control), low N-amended (50 kg ha⁻¹ year⁻¹) and high N-amended (150 kg ha⁻¹ year⁻¹) plots. A total of 1.3 million sequences were processed using qiime. Although Acidobacteria represented the most abundant phylum based on the number of sequences, Proteobacteria were the most diverse in terms of operational taxonomic units (OTUs). UniFrac analyses revealed that the bacterial communities differed significantly among soil horizons and treatments. Microsite variability among the five subplots was also evident. Nonmetric multidimensional scaling ordination of normalized OTU data followed by permutational manova further confirmed these observations. Richness indicators and indicator species analyses revealed higher bacterial diversity associated with N amendment. Differences in bacterial diversity and community composition associated with the N treatments were also observed at lower phylogenetic levels. Only 28–35% of the 6 936 total OTUs identified were common to three treatments, while the rest were specific to one treatment or common to two.

bacterial composition, forest soils, indicator species, pyrosequencing, qiime software, unique operational taxonomic units

Introduction

The composition of soil bacterial communities is influenced by numerous factors ranging from their geographic location (Fulthorpe et al., 2008) to site-specific environment. These factors include, but are not limited to, soil type (Roesch et al., 2007), soil pH, (Neufeld & Mohn, 2005; Fierer & Jackson, 2006; Lauber et al., 2009; Rousk et al., 2010), plant diversity and composition (Carney & Matson, 2006; Uroz et al., 2010), mycorrhizal content (Uroz et al., 2007), and land management practices (Wu et al., 2008; Fierer et al., 2011). Over the last decade, studies of forest soils have provided valuable baseline information about bacterial community structure and diversity in both undisturbed and disturbed forest ecosystems (Axelrood et al., 2002; Chow et al., 2002; Hackl et al., 2004; Chan et al., 2006). Additional studies are needed to address the effects of forestland management practices on bacterial community structures.

Nitrogen (N) is a key regulator of ecosystem processes (Nadelhoffer et al., 1984; Nadelhoffer, 2001; Galloway et al., 2004). Whereas low N availability often limits plant growth and forest productivity in a temperate forest, soil C has been implicated as the key limiting factor for microbial metabolism (Aber et al., 1998; Demoling et al., 2008; Allison et al., 2009, 2010). In the northeastern United States, N concentrations in forest soils have increased over time because of emissions from fossil fuel combustion (Galloway et al., 1984). Bobbink et al. (2010) discussed the worldwide threat of environmental N deposition for eight major ecosystems and its influences on vegetation, including changes in plant communities, plant physiology, and the resistance of plants to pathogens and insect pests. Earlier reports on soils from various forests revealed that N fertilization reduces microbial biomass and activity (Arnebrant et al., 1996; Thirukkumaran & Parkinson, 2000; Bowden et al., 2004; Frey et al., 2004; Wallenstein et al., 2006). A range of inorganic N fertilizer treatments was also demonstrated to decrease soil microbial respiration levels in aspen and pine forests as well as in grasslands (Ramirez et al., 2010).

It has been proposed that elevated N in the soil may cause saturation of the N-retention capacity of forest ecosystems (Aber et al., 1989; Aber, 1992). To study plant and soil microbial responses to elevated N in the soil and to mimic future N additions due to anthropogenic factors, experimental plots subjected to chronic N addition were established in 1989 within pine and mixed hardwood stands at the Harvard Forest Long-Term Ecological Research (LTER) site (http://harvardforest.fas.harvard.edu) in Petersham, MA. These plots have been subjected to NH₄NO₃ addition on a regular basis for over 20 years. Nitrogen saturation due to these treatments has negatively impacted the soils and tree growth at this site, presumably by increasing physiological stress, which leads to tree mortality (Aber et al., 1993, 1998; Magill & Aber, 1998, 2000; Bauer et al., 2004; Magill et al., 2004; McDowell et al., 2004). Frey et al. (2004) and Wallenstein et al. (2006) reported major changes in the bulk composition of soil microorganisms in response to long-term N addition; there was also a decrease in the active fungal biomass, with little change in the bacterial biomass. This resulted in negative shifts in fungal: bacterial biomass ratios, microbial biomass carbon (MBC), and selective induced respiration (SIR) rates, hinting at alterations in the microbial community composition (Wallenstein et al., 2006). The decline in fungal biomass was largely attributed to a decrease in ectomycorrhizae forming associations with major hardwood species, suggesting that fungi are more susceptible to long-term N fertilization than bacteria. However, no details are available concerning the diversity of either the bacterial or fungal populations that may have been impacted by N amendment at this site. The structure of bacterial communities, including N-fixing (Tan et al., 2003), nitrifying, (Compton et al., 2004; He et al., 2007), and denitrifying (Enwall et al., 2005) bacteria, has been confirmed to change in response to N fertilization in forested and agricultural soils. In a previous study by our group, a Ca-supplemented watershed at the Hubbard Brook Experimental Forest (NH) exhibited an altered bacterial composition compared to a reference watershed (Sridevi et al., 2012).

The overall goal of the present study was to test the hypotheses: (1) long-term chronic N-amended organic and mineral soils depict different bacterial community composition compared with control soil; (2) Chronic N-amended soils show an abundance of selected taxa involved in N cycling. Using high-throughput 16S rRNA pyrosequencing (454), we investigated the changes in bacterial community structure in the soils of mixed hardwood stand at Harvard Forest in response to two decades of N fertilization. Specifically, we compared (1) the chemistry of the organic and the mineral soil horizons and (2) the bacterial community structure, diversity, and relative abundance associated with each soil horizon and treatment.

Materials and methods

Site description and experimental design

This work was conducted in soils from the mixed hardwood stand included in the chronic N amendment LTER study at Harvard Forest (http://harvardforest.fas.harvard.edu); the site coordinates are 42° 30′ N and 72° 10′ W. The hardwood stands were approximately 95 years old in 2009, and their history is well documented (Aber et al., 1993). The dominant soils at the study site are typic dystrochrepts of glacial origin and stony to sandy loams with a thick organic horizon (Bowden et al., 2004). The vegetation, climate, average rainfall, site topology, and details of the applied N amendments are described in Bowden et al. (2004) and Magill et al. (2004).

Briefly, the study involved three plots that were established as a part of a larger study initiated in 1989 (Aber et al., 1993), two of which were treated with NH₄NO₃ every year during the growing season, whereas the third was kept as a control. Each plot (30× 30 m) was divided into 36 subplots (each measuring 5 × 5 m). The treatments were as follows: untreated control (Con), low N (LN, treated with 50 kg N ha⁻¹ year⁻¹), and high N (HN, treated with 150 kg N ha⁻¹ year⁻¹). In September 2009, soil cores were collected from five randomly selected subplots within each treatment plot using a soil corer (7.5 cm diameter) to a depth of approximately 15 cm. For each soil core, the upper, dark brown organic horizon (Org) was separated from the lower, light brown mineral horizon (Min). A total of 30 soil samples (five cores per plot × two horizons × three treatments) were collected and brought to the laboratory on ice in polyethylene bags. Each sample was given a designation consisting of the treatment name-subplot name-soil horizon. The samples were sieved (2 mm pore size) to remove roots and stones and stored at −20 °C until further use.

Soil chemical analyses

Air-dried soil samples (20–40 g) were sent to the Soil Testing Service Laboratory at the University of Maine, Orono, ME (http://anlab.umesci.maine.edu) for chemical analyses. The total N and C contents were measured by combustion analysis at 1350 °C. Exchangeable P, Al, and base cations were extracted with 1 M NH₄Cl at a ratio of 2 g of organic soil or 5 g of mineral soil to 100 mL of extraction solution (Blume et al., 1990). After shaking for 1 h, the extracts were vacuum filtered (Whatman # 42, Whatman Inc., Clifton, NJ) and analyzed using flame emission spectroscopy for K and Na and plasma emission spectroscopy for Ca, Mg, and Al. The exchangeable acidity was measured in 1 M of KCl extracts using endpoint titration (Blume et al., 1990). The percentage of soil organic matter in the oven-dried samples was determined using loss-on-ignition (LOI) for 12 h at 550 °C with a muffle furnace. The effective cation-exchange capacity (ECEC) was calculated as the sum of the exchangeable base cations (Ca, Mg, K, and Na) plus the exchangeable acidity.

DNA isolation, PCR amplification, and pyrosequencing

Total DNA was isolated from 0.5 g of each soil sample using the PowerSoil^® DNA isolation kit (MO-BIO Laboratories, Carlsbad, CA) following the manufacturer's instructions, with minor modification that included bead beating for only 7 min with vortex adaptor to minimize DNA shearing. The obtained DNA was quantified and examined for purity (A₂₆₀ : A₂₈₀ ratio between 1.6 and 1.8 and A260 : A230 ratio between 2.0 and 2.2) with a NanoDrop spectrophotometer (Thermo-Fisher Scientific, Waltham, MA). In preparation for pyrosequencing, universal primers with 30 different barcodes (10 bp, one for each of the 30 soil samples) were designed to amplify a 433-bp fragment of the hypervariable region (V6-V8) of the bacterial 16S rRNA gene from the soil samples. The primers used were F968 (5′AA CGC GAA GAACCT TAC3′) and R1401-1a (5′CGG TGT GTA CAA GGC CCG GGA ACG3′) as described in Brons & van Elsas (2008). The amplification reactions were conducted in triplicate using Phusion^®Taq Master Mix (New England Biolabs, Ipswich, MA) with 50 ng of template DNA in a final volume of 50 μL. The reactions were performed in a PTC-100^® Programmable Thermal Cycler (MJ Research, Inc., Waltham, MA) with the following conditions: an initial denaturation at 95 °C for five min, followed by 20 cycles of denaturation at 95 °C for 30 s, annealing at 61 °C for 30 s, and extension at 72 °C for 45 s, with a final extension at 72 °C for 10 min. The triplicate reaction products (amplicons) from each soil sample were pooled for sequencing (Margulies et al., 2005). The pooled PCR products were purified using a DNA purification kit (Zymo Research, Irvine, CA) and subjected to further cleaning via the Agencourt^® AMPure^® XP Bead Purification method (Agencourt Bioscience Corporation, Beverly, MA) to remove fragments < 100 bp. The quality of the PCR products was evaluated in an Agilent 2100 Bioanalyzer using the DNA 1000 LabChip (Agilent Technologies, Palo Alto, CA). The 30 bar-coded samples were pooled in equimolar quantities and processed for sequencing in a full picotiter plate (Roche 454 GS-FLX Titanium System) at the University of Illinois, (www.biotec.illinois.edu/centers/Keck/Highthroughput/).

Data processing

The sequencing data were processed using the Quantitative Insights into Microbial Ecology (qiime) toolkit –Version 1.4.0 (Caporaso et al., 2010b) with default settings for most of the steps, except where specified. The pipeline used for data processing is shown in Supporting information, Fig. S1.

In the first step, the sequences were trimmed and assigned to each soil sample based on their barcodes. Using average quality scores and other parameters, the sequences were also trimmed for quality according to published recommendations (Huse et al., 2008). Multiple steps were required to trim the sequences, including removal of sequences that exhibited the following characteristics: < 200 bp, ambiguous base calling, six-base homopolymer runs, lack of primers, primer mismatches, or uncorrectable barcodes. Similar sequences were clustered into operational taxonomic units (OTUs) using de novouclust (Edgar, 2010) set at a 97% identity threshold, which is a generally acceptable convention for defining an OTU (Lauber et al., 2009). Each OTU consists of a group of sequences identified in a taxonomic study, without designation of their terminal taxonomic rank. The most abundant sequence in the cluster for each OTU was selected as a representative sequence for that OTU. The OTU representative sequences were submitted to NCBI Genbank and were assigned accession numbers JQ049082–JQ060105.

The representative OTU sequences were aligned using pynast (Caporaso et al., 2010b) with the Greengenes database (core set aligned December 16, 2010) (DeSantis et al., 2006), which was set at a minimum sequence identity of 75% and a minimum sequence alignment length of 100 bp. Putative chimeric sequences were removed with Chimera Slayer (Haas et al., 2011). A taxonomic designation was assigned to each representative sequence using the ribosome database project (RDP) classifier with a minimum confidence level of 0.80 (Cole et al., 2005). Singleton OTUs were filtered from the OTU table using the script filter_otus_from_otu_table, and the output was used to perform all other downstream steps. The Lanemask file (the Greengenes core set in the qiime workflow) was used to screen out hypervariable regions, and the default settings were followed for the gap fraction and the alignment threshold in the filtering step. Each of the 30 samples was rarefied to the soil sample exhibiting the lowest number of reads (11 000 sequences) for both α diversity (Rarefaction curves) and β diversity (UniFrac) analyses. Rarefaction curves were generated for Chao1, phylogenetic diversity, and observed species. Alpha rarefaction analyses were performed by computing the average richness metric value from five subplots for each of the six soil types (three treatments × two horizons). Unweighted UniFrac distances were computed based on the Jackknifed OTU table to create a UPGMA tree using the Fast Tree algorithm and visualized with the FigTree programme (Price et al., 2009). This dendrogram was then used to analyze differences in bacterial community composition (Lozupone et al., 2006). Weighted UniFrac distances were used to create a 2D PCoA plot. Processed data were visualized using qiime and pc-ord. Sequence and OTU data were also processed using Microsoft Excel.

Unique OTU calculations

Operational taxonomic unit data were processed separately for each soil horizon and for each treatment within a horizon for (1) the mean of five subplots and (2) the number of OTUs that were unique to each treatment (where each unique OTU was counted only once, regardless of how many of the five subplots contained a copy). The number of OTUs that were unique to the organic or mineral soil samples was also calculated by counting each unique OTU only once regardless of how many of the 15 samples contained copies of that OTU (three treatments × five subplots/treatment). Similarly, the number of OTUs that were unique to all 30 soil samples collected from the Harvard Forest soils was determined by counting each unique OTU only once, regardless of how many of the total 30 soil samples contained a copy of the OTU.

Statistical analyses

systat (Version 10.2) software was used to carry out standard statistical tests, including paired t-tests and two-way analysis of variance (anova) on the soil chemistry data. Nonparametric Kruskal–Wallis test was performed on the taxonomic data. Nonmetric dimensional scaling (NMS) analyses were conducted using pc-ord (Version 6.03, MJM Software Design, Gleneden Beach, OR). To normalize the data, digit one was added to all data before log10 transformation. Bray–Curtis distances were used to delineate patterns between the various treatments and soil horizons (Kruskal, 1964; Mather, 1976). The following conditions were used for the NMS analyses: number of axes = 3, maximum number of iterations = 500, stability criterion (the standard deviation in stress over the last 10 iterations) = 10⁻⁶, number of runs with real data = 100 and the number of runs with randomized data = 250. Random numbers were chosen as a source of starting ordinations. The tie handling was done by penalizing unequal ordination distance (Kruskal's secondary approach). The following were chosen as output options: varimax, randomization test, plot stress vs. iterations, and calculate scores for OTUs by weighted averages. Two-dimensional solutions were finally selected for these analyses based on the assessment using a graph of stress as a function of dimensionality (scree plot). A Monte Carlo test was used to test the stress and the strength of the NMS results. Two-way permutational manova was conducted using the Bray–Curtis distances to evaluate the effect of the horizons and the treatments and the interaction between them. A Mantel test was conducted to evaluate the significance of the correlations among the Bray–Curtis distance scores and the soil chemistry parameters. A joint plot overlay (soil chemical parameters) indicates the relative significance of each axis on a variable based on its length and direction. Indicator species analyses were conducted separately for organic and mineral soil samples using the following setup: treatments as the grouping variable, Dufrêne & Legendre's Indicator Species Analysis (quantitative or presence-absence data) and randomization test with 9999 runs. Dufrene & Legendre's (1997) method combines information on the concentration of OTUs abundance in a particular group with its occurrence in that group. It produces indicator values for each OTU for each treatment. These are then tested for statistical significance using a randomization technique. The indicator data reported here are for P ≤ 0.05.

Results

Soil chemistry

The chemistry of the Con-Org (control organic) soil was significantly different from that of the Con-Min (control mineral) soil (Table 1). The former exhibited a lower pH and higher concentrations of all of the measured inorganic ions, acidity, cation-exchange capacity (CEC), % loss-on-ignition (LOI),% total N, and% total C, but similar C/N ratios. Amendment of the soil with low N did not affect any of the measured parameters, with the exception of resulting in a 30–50% reduction in Na ions in LN-Org (low N organic) and LN-Min (low N mineral) horizons (Table 1). However, the HN-Org (High N organic) soil possessed significantly lower K and Zn concentrations compared with the Con-Org soil; little or no change was observed for the HN-Min (High N mineral) soil.

Soil chemistry of the control and the N–amended soils divided into their organic and mineral horizons

	Organic horizon			Mineral horizon
Mean of five subplots	Control	Low N	High N	Control	Low N	High N
Soil pH	4.00 ± 0.10	4.10 ± 0.10	4.00 ± 0.15	4.60 ± 0.02a	4.70 ± 0.07a	4.60 ± 0.09a
% LOI (organic matter)	28.1 ± 4.6	28.7 ± 4.6	27.2 ± 3.8	9.3 ± 0.4a	9.5 ± 0.8a	11.2 ± 0.4a
Total N (%)	0.60 ± 0.08	0.60 ± 0.04	0.60 ± 0.07	0.20 ± 0.01a	0.20 ± 0.02a	0.20 ± 0.01a
Total C (%)	15.1 ± 2.1	15.4 ± 2.3	15.0 ± 1.9	4.30 ± 0.2a	4.20 ± 0.4a	5.00 ± 0.4a
Ca (mg Kg−1)	115 ± 30	110 ± 12	134 ± 37	13.0 ± 2.4a	16.2 ± 2a	12.6 ± 3a
K (mg Kg−1)	223 ± 54	132 ± 5	118 ± 15**	38.7 ± 9b	37.7 ± 4a	36.3 ± 4a
Mg (mg Kg−1)	102.5 ± 29.1	62.3 ± 6.7	53.6 ± 12.6	11.6 ± 1.9a	10.8 ± 1.9b	12.0 ± 1.1a
P (mg Kg−1)	24.6 ± 8.1	10.1 ± 2.0	11.5 ± 2.2	3.30 ± 0.6a	1.60 ± 0.8a	2.60 ± 0.2a
Al (mg Kg−1)	513 ± 28	650 ± 73	595 ± 35	265 ± 29a	224 ± 42a	333 ± 32a
Fe (mg Kg−1)	85.1 ± 7.3	66.5 ± 9.1	70.5 ± 10.9	12.7 ± 1.7a	7.8 ± 4.0a	19.0 ± 5.9a
Mn (mg Kg−1)	26.8 ± 8.1	11.2 ± 2.1	13.5 ± 3.4	5.80 ± 2.5b	3.60 ± 1.9a	2.70 ± 1.7a
Na (mg Kg−1)	18.9 ± 2.3	10.6 ± 1.9**	14.7 ± 2.5	6.00 ± 0.6a	4.10 ± 0.7*,a	7.90 ± 1.1a
Zn (mg Kg−1)	17.80 ± 4.36	9.40 ± 1.26	7.60 ± 1.04**	2.10 ± 0.69a	2.00 ± 0.24a	1.80 ± 0.28a
Acidity (meq 100 g−1)	9.80 ± 1.10	9.30 ± 1.06	11.10 ± 0.98	4.00 ± 0.28a	3.40 ± 0.64a	4.90 ± 0.56a
CEC (meq)	11.9 ± 0.7	10.8 ± 1.2	12.6 ± 1.3	4.30 ± 0.3a	3.70 ± 0.7a	5.10 ± 0.6a
C : N Ratio	24.7 ± 0.6	26.3 ± 2.3	23.7 ± 1.9	21.0 ± 1.6	20.8 ± 1.6	21.1 ± 1.6

The data displayed are the mean ± SE of five subplots. The superscripts ^aP ≤ 0.05 and ^bP ≤ 0.1 denote significant differences between the Org and Min soil samples irrespective of treatment. *P ≤ 0.1 and **P ≤ 0.05 denote significant differences between the Con and N-amended samples.

Pyrosequencing data analysis

A total of 1 273 206 sequences were obtained from the 30 samples; the sequence length ranged from 200 to 521 bp, with an average read length of 451 bp. Approximately 15% of these sequences were removed by quality filtering. The remaining sequences were distributed between forward primer (FP) and reverse primer (RP) sequences (approximately 42% each). The FP sequence data, for which the mean number of sequences per sample (subplot) was 17 578 ± 658, were used for the detailed analyses presented here. Clustering (by uclust) of all of the FP sequences (527 340) resulted in 13 060 unique OTUs (Table S1). The Org soil samples were represented by 276 693 sequences and the Min soil samples by 250 647 sequences. Greater than 99% of the representative OTU sequences could be aligned with pynast. Chimera Slayer detected a total of 2032 chimeric OTUs from all of the 30 samples combined (15% of the 13 060 total OTUs) which left 11 028 OTUs for downstream processing. Further analysis allowed us to assign 83.7% of the sequences from the Org soil samples and 77.7% of the sequences from the Min soil samples to various bacterial phyla. The remaining sequences (15–20%) could not be classified beyond the bacterial domain and were categorized as ‘unclassified bacteria’. A small number (< 1%) of the OTUs remained unclassified under the bacteria division. After removing 4092 singletons from 11 028 OTUs, the remaining 6936 OTUs were processed further. Processing of the RP data showed only minor differences between the two sets of results (< 10% – details not shown).

Diversity analyses

Clustering of the unweighted UniFrac values for all 30 samples via UPGMA (Fig. 1) revealed that the bacterial composition of Org and Min soil samples was similar to each other by about 85% (UniFrac, P test significance analyses, P < 0.001). The dendrogram also showed that for most plots except for Con-Org, one subplot among five was an outlier (Fig. 1). These analyses did reveal some overlaps in the community composition of three treatment plots within each horizon. Two-dimensional PCoA plots using weighted UniFrac distances for all 30 samples depicted that the principal component 1 explained maximum variation of 15.22%, while principal component 2 explained 6.62% (Fig. 2). Organic soil samples were spatially separated from mineral samples. The figure also showed that soils from subplots for each treatment plot grouped together exhibiting similarity. Control samples for each type of soil horizon were distinctly separated from LN and HN samples.

Tree generated with the pairwise un-weighted UniFrac distances between the samples showing a sample on each node. Each tip represents one sample of a total of 30 (three treatments × two soil horizons × five subplots per treatment) and is denoted by the treatment-subplot-soil horizon; all five subplots of each soil type are represented by a single color. The colored nodes on the left indicate P values and significance level: red = < 0.001, highly significant; yellow = (0.001–0.01) significant; green = (0.01–0.05) marginally significant; and blue (0.05–0.1) suggestive.

Two-dimensional PCoA plot generated using pairwise weighted UniFrac distances from the 30 samples.

Nonmetric multidimensional scaling (NMS; the final stress for the two-dimensional solution was 12.97 with 87 iterations, and the P value for stability was 0.00000) followed by permutational manova (with ordination scores, that is, the Bray–Curtis distance, obtained from the normalized total OTU data) predicted that the bacterial community composition of the N-amended samples was significantly different from the control samples within the corresponding horizons (Fig. S2). In case of mineral soil, there was an overlap between LN-Min and HN-Min bacterial communities. Samples from all five subplots within a treatment plot clustered together and displayed stronger similarities among their OTUs compared to the subplots within other treatment plots. This analysis also revealed significant differences (P < 0.0001) between the Org and Min soil horizons, irrespective of the treatment. This result could also be inferred from the length of the Bray–Curtis distance between the centroids of all of the Org vs. Min subplot cloud points (Fig. S2). The permutational manova results revealed that there was no significant interaction between the horizon and the treatment. Whereas the horizon explained approximately 31% (P ≤ 0.0001) of the variation among the samples, the treatment accounted for < 10.4% (P ≤ 0.002) of the variation.

The Mantel test revealed that there were strong positive correlations of P,K, and acidity with the ordination scores in the Org horizon (Table 2). In the Min horizon soils, only Zn exhibited a significant positive correlation (Table 2). A strong positive correlation with acidity indicated it strongly influenced the bacterial communities in the Org soil samples (Table 2).

Relationship between the soil chemistry and Bray-Curtis (Sorenson) distance measures of the log10-transformed total OTU data (Mantel test)

Soil chemistry	Organic horizon		Mineral horizon
Mantel r	P value	Mantel r	P value
pH	−0.102	0.320	0.202	0.267
Ca	−0.038	0.720	0.176	0.280
P	0.235	0.033**	0.020	0.912
Mn	0.151	0.205	−0.093	0.467
K	0.201	0.091*	0.086	0.622
Zn	0.178	0.098	0.320	0.071*
Na	0.222	0.047	−0.034	0.860
Mg	0.122	0.292	0.050	0.764
Acidity	0.337	0.001**	−0.164	0.291
Al	0.155	0.176	0.029	0.843

All 15 samples from each horizon were pooled for these analyses. The iterations were set to 1000. *P ≤ 0.01 and **P ≤ 0.05 indicate significant correlations.

Richness indices, such as PD, Chao, and Observed species, revealed significantly higher values in the N-amended soils when compared to control soils (P ≤ 0.05) that indicate greater richness of bacterial taxa in N-amended soils (Table 3). Alpha rarefaction curves also showed increasing richness with an increase in the number of sequences analyzed for the N-amended vs. the control samples implying higher richness in N-amended soils (Table 3, Fig. S3).

Effect of N treatment on OTU diversity and richness, as indicated by phylogenetic diversity (PD) estimates, Chao1 values, and the observed species

Treatment	PD estimate mean ± error	Chao1 mean ± error	Observed species mean ± error
Pooled organic and mineral horizon
Con	93 ± 8	2025 ± 249	1294 ± 132
LN	101 ± 7**	2330 ± 170**	1462 ± 90**
HN	102 ± 8**	2320 ± 271**	1482 ± 150**
Organic
Con-Org	90 ± 8	1974 ± 244	1278 ± 121
LN-Org	96 ± 6*	2207 ± 150**	1415 ± 87**
HN-Org	97 ± 7*	2217 ± 261**	1401 ± 145**
Con-Min	95 ± 8	2076 ± 243	1310 ± 141
LN-Min	106 ± 2**	2452 ± 75**	1508 ± 64**
HN-Min	107 ± 5**	2423 ± 239**	1562 ± 106**

The data depicted are the mean of five subplots for each treatment in the organic and/or mineral horizons ± error at 97% confidence intervals. These estimates were calculated in qiime after normalization (rarefaction) of the data to 11 000 sequences per sample. *P ≤ 0.1 and **P ≤ 0.05 denote significant differences between the Con and N-amended samples. A two-sampled t test was conducted in qiime to evaluate differences among the samples.

Unique OTUs

The total number of OTUs that were unique to each treatment (each OTU was counted once, regardless of how many of the subplots contained a copy) was higher than the mean OTU number (Table 4). The sequence abundance varied among the OTUs within each soil horizon: 2% of the OTUs unique to each soil horizon contained 48–52% of the sequences (high abundance OTUs); an additional 18% of the OTUs accounted for 37–41% of the sequences (intermediate abundance OTUs); and > 80% of the OTUs could be classified as low abundance OTUs and these accounted for only 10–14% of the total sequences (Fig. 3). This group of low abundance OTUs was the most diverse in terms of the identifiable bacterial taxa across various phyla (data not shown).

Relationship between the OTUs unique to each treatment and the percent of the total sequences within a treatment for the (a) organic and (b) mineral horizons. Each unique OTU was counted only once, regardless of how many of the five subplots contained a copy.

The number of mean OTUs and the number of OTUs unique to each treatment for the organic and mineral soil horizons

	OTU numbers in the organic horizon
Mean ± SE	Unique to	Mean ± SE	Unique to	Mean ± SE	Unique to
Phylum	Con-Org	Con-Org	LN-Org	LN-Org	HN-Org	HN-Org
Acidobacteria	290 ± 16	535	357 ± 21*	618	321 ± 18	604
Actinobacteria	126 ± 9	236	129 ± 6	226	139 ± 2	249
Bacteroidetes	7 ± 2	26	6 ± 1	21	4 ± 1	16
Chlamydiae	95 ± 7	282	168 ± 25**	451	129 ± 18**	384
Chloroflexi	0 ± 0	2	1 ± 0	3	1 ± 0	3
Firmicutes	10 ± 3	40	10 ± 2	31	14 ± 7	50
Gemmatimonadetes	6 ± 1	15	7 ± 1	16	6 ± 1	11
Nitrospira	0 ± 0	1	1 ± 0*	1	1 ± 0*	1
Proteobacteria	435 ± 32	1066	582 ± 48*	1362	513 ± 15	1223
TM7	16 ± 4	51	20 ± 3	56	25 ± 1*	68
Verrucomicrobia	118 ± 7	204	130 ± 6	219	131 ± 5	1097
Un-Bacteria	357 ± 28	863	464 ± 41*	1088	472 ± 36**	225
	OTU numbers in the mineral horizon
Mean ± SE	Unique to	Mean ± SE	Unique to	Mean ± SE	Unique to
Phylum	Con-Min	Con-Min	LN-Min	LN-Min	HN-Min	HN-Min
Acidobacteria	285 ± 13	517	326 ± 16	572	322 ± 14*	571
Actinobacteria	99 ± 12b	200	114 ± 13	230	118 ± 6	222
Bacteroidetes	4 ± 1	18	4 ± 1	14	5 ± 1	14
Chlamydiae	114 ± 17	357	170 ± 8**	470	143 ± 22	436
Chloroflexi	3 ± 0a	4	3 ± 0	4	3 ± 1	4
Firmicutes	29 ± 3a	77	30 ± 8	77	20 ± 4	55
Gemmatimonadetes	6 ± 0	10	8 ± 0**	14	7 ± 1	10
Nitrospira	1 ± 0a	2	1 ± 0	3	1 ± 0	2
Proteobacteria	386 ± 42	1042	495 ± 43	1269	451 ± 34	1131
TM7	11 ± 2	39	13 ± 2	47	10 ± 1	40
Verrucomicrobia	109 ± 6	199	126 ± 8	239	120 ± 6*	226
Un-Bacteria	468 ± 30	1073	583 ± 26**	1294	585 ± 30**	1262

Each unique OTU was counted only once for the organic or mineral horizon of each treatment, regardless of how many of the five treatment subplots contained a copy of the OTU. Un-bacteria denotes unclassified bacteria. ^aP ≤ 0.05 and ^bP ≤ 0.1 denote significant differences between the Org and Min soil samples, irrespective of the treatment. **P ≤ 0.05 and *P ≤ 0.1 denote significant differences between the Con and N-amended samples.

Con-Org and Con-Min contained 3324 and 3543 unique OTUs, respectively, with approximately 45% of these OTUs being shared between the two horizons (Fig. 4a). Therefore, only 4729 unique OTUs were found in the control soil samples (Org plus Min) across both horizons.

Venn diagrams showing the distribution of the operational taxonomic units (OTUs) that were unique or common to the (a) Control Organic and Control Mineral soil horizons; (b) Organic horizon for Con, LN, and HN treatments; and (c) Mineral horizon for the three treatments.

The number of OTUs in the LN (4096) and HN (3937) samples was greater than in the control (3324) samples (Fig. 4b). The number of unique OTUs within the Org soils was 5711 of a total of 11 357 identified from all 15 samples. For this calculation, each OTU was counted once, regardless of how many of the 15 samples contained a copy. Similarly, in the mineral soil horizon, the LN (4238) and HN (3979) treatments exhibited higher numbers of unique OTUs compared with the control samples (3543) (Fig. 4c). The number of unique OTUs within the Min soils was 5836 of a total of 11 760 identified from all 15 samples.

Further examination of the organic and mineral horizon OTU data revealed that among the 5711 total OTUs unique to the 15 organic samples, 9, 13, and 14% were found only in the Con-Org, LN-Org, and HN-Org samples, respectively (Fig. 4b); among the 5836 OTUs unique to the Min horizon, 9, 14 and 12% were found only in the Con-Min, LN-Min, and HN-Min samples, respectively (Fig. 4c). From 28 to 35% of the OTUs were common to all of the treatments within each horizon, whereas only approximately 10% were common to the two N-amended soils within each horizon. Overall, a greater number of OTUs were common to the N-amended soils (13–14%) than to either the LN or the HN samples and the control (5–9%) samples for each horizon (Figs. 4b and 4c). The OTUs that were unique to each treatment were further processed at the phylum level, and after the OTUs common among the three treatments were removed, the net result was the number of OTUs that were specific to each treatment (Fig. 5a and b). These data revealed that the number of treatment-specific OTUs was higher for the N treatments. Thus, the diversity and the composition of the N-amended soil were higher compared to the control soil. More details regarding the unique and shared OTU distributions among the treatments at the phylum level are shown in Table S2.

The effect of N treatment on the OTU numbers per phylum for the organic (a) and mineral (b) soil horizons. Each bar represents a unique set of OTUs that are only present in a particular sample (e.g. the OTUs present in Con-Org are absent from the N-amended LN-Org or HN-Org samples), thus indicating the change in diversity and composition with N amendment. OTUs shared among more than one treatment are absent from this figure.

Bacterial community composition of untreated (control) soils of Harvard Forest

Sequence tables were created by counting only those genera that were represented by 10 or more mean sequences in at least one soil type of six. The phylum Acidobacteria was represented by the greatest number of sequences (43–48%) in both horizons of the control soils at this site (Table 5). This group was followed in abundance by the Proteobacteria, Actinobacteria, and Verrucomicrobia, which ranged from 4 to 16% in terms of relative sequence abundance (Table 5). The phyla Chlamydiae,Firmicutes, Gemmatimonadetes, TM7, Bacteroidetes, Chloroflexi, and Nitrospira were each contributed ≤ 1–2% of the total sequences. Although the phylum Acidobacteria was the most abundant in terms of sequence number (Table 5), the phylum Proteobacteria was the most diverse in terms of OTUs (Table 4). Whereas there were many OTUs corresponding to the Actinobacteria in the Con-Org samples, the Con-Min samples showed a high number of OTUs representing the Chloroflexi,Firmicutes, and Nitrospira (Table 4).

Effect of N treatment on the relative abundance of sequences within bacterial phyla

Phylum	Organic horizon			Mineral horizon
Con	Low N	High N	Con	Low N	High N
Acidobacteria	7369 ± 385 (48%)	9069 ± 1212	8427 ± 521	6757 ± 254 (43.5%)	8020 ± 584	6164 ± 740
Actinobacteria	1453 ± 179 (9.5%)	1517 ± 141	1636 ± 166	696 ± 222a (10.3%)	950 ± 313	671 ± 73
Bacteroidetes	10 ± 3 (0.07%)	8 ± 2	6 ± 1	9 ± 5 (0.06%)	4 ± 1	8 ± 3
Chlamydiae	321 ± 47 (2.1%)	531 ± 102*	451 ± 61	344 ± 73 (2.2%)	531 ± 64*	393 ± 92
Chloroflexi	1 ± 0 (0.01%)	1 ± 1	2 ± 2	10 ± 2a (0.06%)	15 ± 7	11 ± 2
Firmicutes	19 ± 6 (0.1%)	11 ± 2	33 ± 18	92 ± 28a (0.5%)	110 ± 61	47 ± 13
Gemmatimonadetes	31 ± 6 (0.2%)	36 ± 5	32 ± 5	62 ± 21 (0.4%)	78 ± 45	65 ± 16
Nitrospira	0 ± 0 (0%)	2 ± 1	1 ± 1	6 ± 1a (0.04%)	42 ± 15**	16 ± 10
Proteobacteria	2355 ± 355 (15.4%)	3544 ± 335**	3339 ± 278*	2270 ± 679 (14.6%)	2646 ± 358	2274 ± 285
TM7	31 ± 12 (0.2%)	32 ± 7	64 ± 8**	21 ± 4 (0.1%)	18 ± 3	14 ± 2
Verrucomicrobia	1286 ± 95 (8.4%)	1418 ± 175	1827 ± 196**	2105 ± 258a (13.5%)	2153 ± 146	1568 ± 79*
Un-Bacteria	2325 ± 217 (15.2%)	3016 ± 440	3411 ± 246**	3140 ± 285a (20.2%)	4011 ± 693	3761 ± 160

The data depicted are the mean ± SE of five subplots. Un-Bacteria denotes unclassified bacteria. The values in parentheses represent a percent of a given category within the total number of bacteria; percent data are presented only for the control treatments. This table was created by counting only those genera that were represented by 10 or more mean sequences in at least one soil-type out of six. ^aP ≤ 0.05 and ^bP ≤ 0.1 denote significant differences between the Org and Min soil samples, irrespective of the treatment. *P ≤ 0.1 and **P ≤ 0.05 denote significant differences between the Con and N-amended samples.

In terms of relative sequence abundance, no significant differences were observed between the Con-Org and the Con-Min samples for Acidobacteria, Bacteroidetes,Chlamydiae, TM7, and Proteobacteria (Table 5). Other noteworthy differences in terms of the relative sequence abundance at a lower phylogenetic level are shown in Table S3.

Bacterial communities differed among the treatments

At the phylum level, significantly more Chlamydiae and Proteobacteria sequences were observed in the LN-Org samples compared to the control samples (Table 6), whereas for TM7 and unclassified bacteria sequences were higher in the HN-Org samples compared to the Con-Org. No differences with respect to sequence numbers were observed for the mineral soils at the phylum level. In terms of OTUs, LN-Org exhibited greater numbers corresponding to Acidobacteria,Chlamydiae,Nitrospira,Proteobacteria, and unclassified bacteria (Table 4). Among the mineral samples, the phyla Chlamydiae and Gemmatimonadetes displayed higher numbers of OTUs in the LN-Min samples. Additionally, Acidobacteria and Verrucomicrobia OTUs were more abundant in the HN-Min samples, whereas unclassified bacteria were more abundant in both types of N-amended Min soil samples (Table 4). To assess the differences in bacterial community structure between the control and the N-amended samples (i.e. the mean OTU and sequence numbers) at lower phylogenetic levels, analyses were performed at genus level (Table S3, Fig. S4). The OTUs and sequences corresponding to unclassified bacteria were more abundant in the N-amended samples vs. the controls from both horizons (Tables 4 and 5).

Unique indicator OTUs present in each soil type and treatment. OTUs that were unique to a plot at P ≤ 0.05 are listed below in this table. For classification of each indicator OTUs at genus level refer to Tables S3, S4 and S5

Soil Horizon	Unique indicator Species present in each soil horizon
Organic	176
Mineral	186
	Unique indicator species comparison among three types of treated soils within each horizon
Soil Horizon	Con	LN	HN
Organic	27	82	77
Mineral	21	40	59

Indicator OTUs

Comparison of organic and mineral soil

Indicator species analyses revealed that a total of 176 and 186 OTUs were significant indicators of organic and mineral horizon, respectively (P ≤ 0.05, Table 6). In phylum Acidobacteria, specific indicator OTUs belonging to subdivisions Gp1, 2, 3 and 13 could be found in each horizon but those of Gp 4, 5, 6, 7 &10 were present only in mineral soils (Table S4). Phylum Actinobacteria was represented by indicator OTUs of genus Sporichthya,Actinospica and Conexibacter in organic soils. In mineral soils it was represented mostly by OTUs of unclassified genera of families Streptomycetaceae and Rubrobacteraceae. While indicator OTUs of genus Gemmatimonas (phylum Gemmatimonadetes) were present in each horizon, presence of OTUs belonging to Parachlamydia (Chlamydiae) was indicative of mineral soils. Specific OTUs of Chloroflexi,Firmicutes, and Nitrospira were indicative of mineral soils. In Proteobacteria while some indicator OTUs of classes α, δ & γ were present in each horizon, OTUs of β-Proteobacteria were indicators only of mineral soils. In phylum Verrucomicrobia the OTUs of genera Opitutus and Verrucomicrobium were linked with organic soils while that for Subdivision 3 and Xiphinematobacteriaceae were present in both horizons. Twice as many unclassified indicator OTUs were found in mineral horizon compared to organic.

Comparison of three treatments within each soil horizon

In organic soils, in phylum Acidobacteria specific indicator OTUs of subdivision Gp1 were present in each of the three organic soil plots whereas those of Gp2 were present only in N-amended soils (Table S5). HN-Org, on the other hand, had indicator OTUs of Gp13 and Gp6 subdivisions. Genus Opitutus (Verrucomicrobia) was represented by indicator species in Con-Org and HN-Org soils (Table S4). Specific OTUs of genera Actinospica (Actinobacteria), Mycobacterium (Actinobacteria), Parachlamydia (Chlamydiae), Rhabdochlamydia (Chlamydiae), and Bacillus (Firmicutes) were indicators of LN-Org soils. Indicator OTUs of genus Aquicella (γ-Proteobacteria) were found in LN-Org and HN-Org soils.

In mineral soils, specific indicator OTUs for subdivisions, Gp1 and 3 (Acidobacteria) were present in each of the three treatments; indicator OTUs of Gp 6 were present only in LN-Min soils and for Gp5,10, and 13 were present only in HN-Min soils (Table S6). Genus Opitutus (Verrucomicrobia) was represented by specific indicator species in all three soils. Indicator OTUs representing the genus Sporichthya (Actinobacteria) were present in LN and HN-amended mineral soils. Specific OTUs of Nitrospira were indicative of LN-Min soils (Table S6). For more details on indicator species data see the Tables S4 and S5.

Discussion

The present study tested whether N-amended forest soils had altered bacterial community composition compared to control soil and that these alterations were specific to each soil horizon. Data for zero time are not available for this study because high-throughput analytical tools for analysis of microbial populations did not exist in 1989. In addition, none of the groups working at present at Harvard Forest have frozen soils for time zero. This has also been the case for other studies that used similar approaches to analyze microbial populations in different environments (Campbell et al., 2010; Shen et al., 2010; Fierer et al., 2011; Sridevi et al., 2012).

UPGMA tree generated based on the unweighted UniFrac distances revealed significant differences between two soil horizons for each of the three treatment plots. This analysis also revealed that there were differences among the five subplots within each treatment plot, which can be ascribed to microsite-dependent heterogeneity in soil chemistry. The higher number of unique OTUs compared to the mean number of OTUs found for each treatment (Table 4) supports this argument. These results are consistent with those of Sridevi et al. (2012), who also observed microsite variation in soil samples from the Hubbard Brook Experimental Forest, NH. However, the variability among the subplots was much lower than the variability between the treatment plots, as demonstrated by the weighted PCoA and weighted NMS analyses.

Campbell et al. (2011) reported that in coastal bacterial populations, a small number of OTUs contained the majority of sequences, whereas a limited number of sequences were present in most OTUs. Based on the comparison of rRNA activity with the respective OTU frequencies, these authors further concluded that the rare group of OTUs represented more active bacteria compared to abundant OTUs. In the present study, approximately 2% of the OTUs contained 50% of the total sequences in each sample, and 10% of the total sequences generated 80% of the OTUs. However, no attempt was made in this study to relate the activity of these OTUs to their relative frequencies.

The majority (70–80%) of the sequences included in these data fell into four known bacterial phyla, viz., the Acidobacteria,Proteobacteria,Actinobacteria, and Verrucomicrobia, all of which are considered to be universally present in soils (Hugenholtz et al., 1998; Janssen, 2006). The phylum Acidobacteria was the most abundant in terms of sequence numbers, followed by the Proteobacteria, which is in agreement with previous reports on forest soils across different climatic conditions (Janssen, 2006; Jones et al., 2009; Nacke et al., 2011). It has been suggested that soils sharing similar physicochemical characteristics harbor similar bacterial communities, despite geographic distances (Fierer & Jackson, 2006). According to Smit et al. (2001), the Acidobacteria to Proteobacteria ratio of > 1.0 may indicate low nutrient conditions in the soils of the Harvard Forest. Among the Acidobacteria, 10 subgroups (1, 2, 3, 4, 5, 6, 7, 8, 10 and 13) were detected out of a total 26 groups identified thus far (Hugenholtz et al., 1998; Zimmermann et al., 2005; Barns et al., 2007). Among the Proteobacteria present, the sequences from the γ-proteobacteria class were more abundant than those of the other three classes.

The bacterial communities in the organic and the mineral soil horizons are different

Whereas the organic and the mineral horizons of the control (untreated) soil samples exhibited similar total numbers of OTUs, the bacterial community structure within the two horizons was significantly different as demonstrated by the weighted NMS and weighted PCoA analyses. Only 25–30% of the total OTUs were specific to each horizon which indicates that the two horizons harbored substantially different bacterial taxa despite vertical separation of only a few centimeters. This fact was also validated by indicator species analyses. Differences in the soil chemistry of the two horizons are likely to be the major contributors to this variation. Axelrood et al. (2002) reported major differences in the bacterial composition of organic and mineral soils from forests in British Columbia, where Actinobacteria dominated the mineral soil in terms of sequence numbers, and the β-proteobacteria class dominated organic soil. However, in the present study, the phylum Acidobacteria was the most abundant in both horizons in terms of sequence numbers. Differences in the aboveground plant communities present at these study sites as well as the different methodologies used (clonal library sequencing vs. pyrosequencing) may account for these divergent results. At Harvard Forest, whereas the diversity (numbers of OTUs) of the Con-Org samples was greater than that of the Con-Min samples for Actinobacteria, the reverse was the case for Firmicutes (Table 4). Twice as many unclassified bacteria were identified as indicator OTUs in mineral soil compared to organic validating prior reports that bacterial taxa of subsurface soils are not as well annotated as the surface soils.

Verrucomicrobia are found more often in deeper soil layers because they show a metabolic preference for low nutrient availability or low C concentrations (da Rocha et al., 2010). Fierer et al. (2003) reported a decrease in soil organic C content with increasing depth. The existence of a greater number of Verrucomicrobia sequences in the Con-Min soil samples compared to the Con-Org samples at Harvard Forest is in accord with the above findings (Table 1).

The bacterial communities are different in the control and the N-amended soils

The observation that only 28–35% of the unique OTUs were common to all three treatments within each soil horizon points to major rearrangements in the bacterial community structure associated with N amendment. Although each of the soil types retained its own unique bacterial community, the LN and HN samples shared a greater number of OTUs than either sample shared with the control samples. The disappearance of some OTUs and the appearance of a larger number of new ones in N-amended soils validates the above point of genera rearrangements in these soils. Indicator species (OTU) analyses validated these unique OTU observations by showing the presence of indicator OTUs that were significantly associated with a particular treatment. This observation also emphasizes the important role of N amendment in shaping the bacterial populations in forest soils. Similarly, at Hubbard Brook Experimental Forest, NH, Ca-supplemented soils exhibited a significantly different bacterial composition (20%) compared to the composition of reference watershed soils (Sridevi et al., 2012).

Although the specific physiological functions of the Acidobacteria in the soil remain unknown, their abundance suggests that this group must play an important role in soil ecosystem functioning. Acidobacteria tend to favor soils with low C availability (Fierer et al., 2007). The dissolved organic C content in the Harvard Forest soils was low regardless of the N amendments applied (Magill & Aber, 2000), which might be responsible for the abundance of this group in soils under all three treatments. The long-term N-fertilization regime also had no effect on the predominant Gp1 and Gp2 subdivisions within the Acidobacteria in terms of their relative abundance and diversity at this site. The proportions of the other subdivisions varied with the applied treatments.

Whereas the Acidobacteria contributed the highest number of sequences in the N-amended samples, the Proteobacteria displayed the greatest number of mean and unique OTUs which implies that the latter group of organisms may adapt better to changes in the soil chemistry. The classes α- and β-Proteobacteria, whose members were most abundant in the N-amended soils at Harvard Forest, contain well-known N-fixing and N-transforming genera (e.g. Bradyrhizobium, Burkholderia, Mesorhizobium, and Magnetospirillum) as well as methane-oxidizing genera (e.g. Methylocella) (Dedysh et al., 1998). In fact, the genus Burkholderia was detected as an indicator species in the HN-Min soil. Unlike previous studies that reported either positive or negative impacts of N fertilizers on the abundance of methanotrophs (Bodelier et al., 2000; Mohanty et al., 2006), the present data did not reveal any change in the numbers of these bacteria. A preliminary study of changes in ammonia monooxygenase (amoA) gene restriction fragment length polymorphism (RFLP) patterns in response to N amendment at Harvard Forest (Compton et al., 2004) found a higher abundance of Nitrospira compared to other ammonia-oxidizing bacteria in N-amended soils. Quantitative PCR analysis of amoA genes also revealed higher abundance of Nitrospira in N-fertilized soils compared to controls (Wertz et al., 2012). Similarly, in our N-amended soil samples, higher numbers of Nitrospira sequences were observed, which indicates a possible increase in nitrification activity. Genus Nitrospira was detected in the LN-Min as an indicator species.

Nemergut et al. (2008) reported a higher relative abundance of the Bacteroidetes and the Gemmatimonadetes and a lower abundance of the Verrucomicrobia (which include methane oxidizers) (Dunfield et al., 2007) in long-term N-treated soils of the alpine tundra. At Harvard Forest, the relative abundance of sequences and unique OTUs for the phylum Verrucomicrobia was higher in the N-amended soils, whereas there was no major change in abundances for the Bacteroidetes and Gemmatimonadetes. Because the two ecosystems possess such different plant communities, these differing responses to N amendment are understandable. In addition, the functional diversity among the members within any individual phylum is typically large; thus, it cannot be expected that all of the members of a phylum will share common ecological and metabolic characteristics. At Harvard Forest, specific OTUs of genus Opitutus (Verrucomicrobia) were detected as indicators of different soil types. This observation reveals that there were major rearrangements of species within this genus with N amendments.

Soil chemistry and bacterial composition are strongly correlated

Whereas reductions in the soil Ca, Mg, and P contents were observed at Harvard Forest from 1995 to 2000 (Minocha et al., 2000), the data included in the current study (from 2009) showed no change in the exchangeable soil Ca and Mg contents but did reveal a decrease in the K, Na, Zn, and Mn contents. Previous studies that have documented the impact of K content on plant growth and physiology support the conclusion that K plays an important role in forest primary productivity (Tripler et al., 2006). The base cation losses in the soil in the N-amended plots at Harvard Forest have shown to be positively correlated with foliar nutrient imbalances, increases in foliar free putrescine (a general stress indicator in plants), growth declines, and increased mortality of red maple (Aber et al., 1995; Minocha et al., 1997, 2000; Bauer et al., 2004; Magill et al., 2004). Maple mortality in the N-amended plots during the past few years and the increase in litter fall in response to N additions may have contributed to the observed changes in the soil Ca, Mg, and P contents in 2009 and related bacterial community changes. Chemical interactions between plant and bacterial communities have been demonstrated for symbiotic associations (Bright & Bulgheresi, 2010 and references therein; Rooney et al., 2010). The positive correlation between the base cations and the Bray–Curtis (Sorensen) distances for the bacterial communities of both soil horizons indicates that some of the variance in the bacterial communities may be caused by the soil chemistry of the treatment plots. These results further confirmed a strong correlation between the soil microflora and changes in soil acidity due to various environmental and anthropogenic factors and fertilizer treatments. Major rearrangements in bacterial communities in N-amended soils without any change in total C and N is an unexpected finding. The exact sequence of events accompanying or causing the changes in soil chemistry and the plant and microflora as well as the ecological and physiological impacts of these changes following N amendments are, however, not known at present.

Conclusions

At Harvard Forest, chronic N treatment affected the diversity and the composition of bacterial communities within organic and mineral soil horizons. N-amended samples exhibited a greater richness of bacterial taxa compared with control samples. Whereas about 50% of the total sequences corresponded to only 2% of the total OTUs, approximately 10% of the total sequences were highly diverse and contributed up to 80% of the total OTUs. Shifts in the bacterial communities were obvious at the family and genus levels. The changes observed in the bacterial community structure in response to chronic N treatments may be a cumulative outcome of N-driven soil base cation changes, net changes in aboveground plant productivity, and changes in fungal biomass as reported earlier by our and other groups. Future studies must examine the metabolic functions of these bacterial populations following changes in response to chronic N amendment and changes that occur after short-term N amendment.

Acknowledgements

We thank Stephanie Long, Gloria Quigley, and Kenneth Dudzik for help in the collection of the soil samples and for editing of the manuscript. We also acknowledge the help of Feseha Abebe Akele in designing the primer bar codes, William A. Walters for help with QIIME, and Dr. Kevin Smith and Dr. Greg Caporaso for reviewing the manuscript. The Maine Soil Testing Service at the University of Maine is acknowledged for their help in the soil analysis. We thank Serita Frey and Scott Ollinger for providing access to the Chronic Nitrogen Amendment Study at Harvard Forest which has been maintained since 1988 by John Aber and others with funding from the National Science Foundation's Long-term Ecological Research (LTER) Program. This paper is scientific contribution Number 2470 from the New Hampshire Agricultural Experiment Station. The authors have no conflict of interest to declare.

Supporting Information

Additional Supporting Information may be found in the online version of this article:

Fig. S1. A detailed flow chart of the third party tools used within the qiime toolkit.

Fig. S2. NMS ordination calculated from OTU raw data from the 30 soil samples. Each soil type is represented by five replicates and a centroid, which is indicated by a single color and number.

Fig. S3. The rarefaction curves of the control and N-amended soil samples of organic (a, c, and e) and mineral horizons (b, d, and f).

Fig. S4. The fold changes in the numbers of OTUs among various known genera and unclassified family members from the N-amended soils relative to the control soils.

Table S1. Raw data on total number of sequences and OTUs for each of the 30 soil samples before and after filtering singletons.

Table S2. Effects of N treatments on relative abundance of OTUs (% of total) of bacterial phyla those were unique to organic or mineral soil horizons.

Table S3. Effect of N treatment on the relative abundance of sequences within Acidobacteria groups and genera of other phyla for organic and mineral soil horizons.

Table S4. Indicator species unique to each untreated soil horizon.

Table S5. Indicator species unique to each treatment in the organic soil horizon.

Table S6. Indicator species unique to each treatment in the mineral soil horizon.

Please note: Wiley-Blackwell is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

Author notes

Editor: Angela Sessitsch