Nutrient addition dramatically accelerates microbial community succession.

The ecological mechanisms driving community succession are widely debated, particularly for microorganisms. While successional soil microbial communities are known to undergo predictable changes in structure concomitant with shifts in a variety of edaphic properties, the causal mechanisms underlying these patterns are poorly understood. Thus, to specifically isolate how nutrients--important drivers of plant succession--affect soil microbial succession, we established a full factorial nitrogen (N) and phosphorus (P) fertilization plot experiment in recently deglaciated (∼3 years since exposure), unvegetated soils of the Puca Glacier forefield in Southeastern Peru. We evaluated soil properties and examined bacterial community composition in plots before and one year after fertilization. Fertilized soils were then compared to samples from three reference successional transects representing advancing stages of soil development ranging from 5 years to 85 years since exposure. We found that a single application of +NP fertilizer caused the soil bacterial community structure of the three-year old soils to most resemble the 85-year old soils after one year. Despite differences in a variety of soil edaphic properties between fertilizer plots and late successional soils, bacterial community composition of +NP plots converged with late successional communities. Thus, our work suggests a mechanism for microbial succession whereby changes in resource availability drive shifts in community composition, supporting a role for nutrient colimitation in primary succession. These results suggest that nutrients alone, independent of other edaphic factors that change with succession, act as an important control over soil microbial community development, greatly accelerating the rate of succession.

Introduction

Deglaciated forefields have been valuable model systems for developing and tesn class="Chemical">ting theories of successioclass="Chemical">n aclass="Chemical">nd have greatly eclass="Chemical">nhaclass="Chemical">nced our uclass="Chemical">nderstaclass="Chemical">ndiclass="Chemical">ng of the relatioclass="Chemical">nship betweeclass="Chemical">n commuclass="Chemical">nity structure aclass="Chemical">nd fuclass="Chemical">nctioclass="Chemical">n duriclass="Chemical">ng ecosystem developmeclass="Chemical">nt [1]–[3]. Shifts iclass="Chemical">n soil class="Chemical">nutrieclass="Chemical">nt pools, iclass="Chemical">ncludiclass="Chemical">ng iclass="Chemical">ncreases iclass="Chemical">n available class="Chemical">n class="Chemical">nitrogen (N) and phosphorus (P), have been well documented along early primary successional chronosequences [4]–[6] and have been shown to correlate with changes in plant community succession [1], [7], [8]. Recently, studies in such systems have revealed that – like plants – microbial communities also progress through successional stages [9]–[11]. However, the forces that control microbial succession are not well understood.

Some evidence suggests that shifts in nutrient availability may also, in part, drive microbial community succession. For example, in primary successional ecosystems, research has corroborated relationships between natural gradients in soil nutrients and microbial community composition [12], [13]. Such correlations can be difficult to interpret, however, as changes in microbial community composition could be both a cause and consequence of shifts in soil fertility. Furthermore, the mechanisms underlying correlations between standing nutrient pools and microbial communities may be temporally disconnected, in that current soil biogeochemical status may not accurately reflect the historical nutrient conditions that structured the microbial community. Thus, manipulation experiments are essential in evaluan class="Chemical">ting the direct impact of class="Chemical">nutrieclass="Chemical">nts aclass="Chemical">nd their limitatioclass="Chemical">ns oclass="Chemical">n microbial commuclass="Chemical">nities. Iclass="Chemical">ndeed, fertilizer treatmeclass="Chemical">nts are kclass="Chemical">nowclass="Chemical">n to elicit chaclass="Chemical">nges iclass="Chemical">n soil microbial commuclass="Chemical">nity structure aclass="Chemical">nd fuclass="Chemical">nctioclass="Chemical">n iclass="Chemical">n more developed ecosystems [14], [15] suggesclass="Chemical">n class="Chemical">ting that nutrient availability may also be important in controlling successional changes in microbial community composition.

Yet, it would be surprising if nutrients alone drove microbial community succession for several reasons. First, other edaphic properties also undergo concomitant shifts with microbial community structure and function during succession, some of which are known to more strongly correlate with microbial community structure than nutrient pools in developed soils. For example, organic n class="Chemical">carbon (C) pools aclass="Chemical">nd pH, which typically show dramatic chaclass="Chemical">nges across primary successioclass="Chemical">nal chroclass="Chemical">nosequeclass="Chemical">nces [2], are key determiclass="Chemical">naclass="Chemical">nts of soil microbial commuclass="Chemical">nity compositioclass="Chemical">n at regioclass="Chemical">nal to global scales [16]–[19]. Secoclass="Chemical">nd, soil microbial commuclass="Chemical">nity structure caclass="Chemical">n correlate with placlass="Chemical">nt commuclass="Chemical">nity compositioclass="Chemical">n [20], [21], which caclass="Chemical">n show stroclass="Chemical">ng spatial gradieclass="Chemical">nts iclass="Chemical">n early successioclass="Chemical">n [3]. Third, stochastic processes caclass="Chemical">n be key iclass="Chemical">n shapiclass="Chemical">ng early successioclass="Chemical">nal commuclass="Chemical">nities where the importaclass="Chemical">nce of dispersal eveclass="Chemical">nts may be acceclass="Chemical">ntuated, [22]–[24] aclass="Chemical">nd arrival order may iclass="Chemical">nflueclass="Chemical">nce assembly through priority effects [25]. Giveclass="Chemical">n the large fuclass="Chemical">nctioclass="Chemical">nal aclass="Chemical">nd phylogeclass="Chemical">netic diversity of microbial commuclass="Chemical">nities, it is possible that successioclass="Chemical">n is iclass="Chemical">nflueclass="Chemical">nced by a diverse combiclass="Chemical">natioclass="Chemical">n of such factors [26].

Thus, the extent to which nutrients themselves influence microbial community assembly outside of the myriad of factors that change over succession is unknown. To specifically isolate the effects of nutrients, we performed a full factorial n class="Chemical">N×P fertilizatioclass="Chemical">n experimeclass="Chemical">nt iclass="Chemical">n soils that had beeclass="Chemical">n exposed for ∼3 years iclass="Chemical">n the forefield of the Puca Glacier iclass="Chemical">n Southeasterclass="Chemical">n Peru. We aclass="Chemical">nalyzed soil bacterial commuclass="Chemical">nities before aclass="Chemical">nd oclass="Chemical">ne year followiclass="Chemical">ng class="Chemical">nutrieclass="Chemical">nt additioclass="Chemical">ns aclass="Chemical">nd compared them with soils sampled from three differeclass="Chemical">nt locatioclass="Chemical">ns over aclass="Chemical">n 85-year sectioclass="Chemical">n of the Puca Glacier chroclass="Chemical">nosequeclass="Chemical">nce. The Puca Glacier soils coclass="Chemical">nstitute aclass="Chemical">n autotrophic successioclass="Chemical">nal sequeclass="Chemical">nce [27], aclass="Chemical">nd both photosyclass="Chemical">nthesis aclass="Chemical">nd respiratioclass="Chemical">n respoclass="Chemical">nd stroclass="Chemical">ngly to P additioclass="Chemical">ns iclass="Chemical">n microcosms [28], [29]. class="Chemical">n class="Chemical">Nitrogen appears to be limiting in this system as well and N-fixation rates in 4 year old unvegetated soils are comparable to rates measured in developed soil crusts [30]. Thus, given work that demonstrates relationships between nutrients and microbial community composition, we hypothesized that fertilizer additions to early successional soils would drive communities to be compositionally different than unfertilized (control) soils. However, due to the potential influence of other edaphic (e.g. pH, organic C, soil moisture) and stochastic factors on microbial succession, we hypothesized that fertilized communities would be unique from communities found along the natural chronosequence.

Materials and Methods

Study site description, fertilization, and sampling

The study site is located in the forelands of the Puca Glacier in the n class="Chemical">Cordillera Vilcanota of Peru (13°46′24″S, 71°04′17″W, ∼5,000 m.a.s.l.). class="Chemical">n class="Chemical">No specific permits were required for our field studies and our work did not involve endangered or protected species. Mean annual precipitation is roughly 100 cm and mean annual temperature is ∼5°C. Moraine rocks at this site have high quartz and calcite mineral content. Further details of this site can be found in previous work [9], [30] and soil characteristics are presented in Table S2.

We established permanent plots (1 m2) near the terminus of the glacier, in soils that had been deglaciated for approximately 3 years at the time of initial sampling. Corners were marked with long nails (approximately 15 cm shank length) to guide resampling. Sampling occurred in August 2010 (pre-treatment) and August 2011 (post-treatment). All of the plots were unvegetated and no n class="Species">mosses aclass="Chemical">nd licheclass="Chemical">ns were preseclass="Chemical">nt at the time of establishmeclass="Chemical">nt. Each of the 16 plots was raclass="Chemical">ndomly choseclass="Chemical">n to receive oclass="Chemical">ne of three class="Chemical">nutrieclass="Chemical">nt ameclass="Chemical">ndmeclass="Chemical">nts (class="Chemical">n class="Chemical">nitrogen addition (+N), phosphorus addition (+P), the combination of the two (+NP)) or to serve as controls, resulting in a total of four plots per treatment and four control plots.

Pre-weighed amounts of fertilizer were dissolved in glacier-melt stream n class="Chemical">water aclass="Chemical">nd fertilizer solutioclass="Chemical">ns were applied with haclass="Chemical">ndheld sprayers. Each sprayer was desigclass="Chemical">nated for a particular treatmeclass="Chemical">nt to avoid cross coclass="Chemical">ntamiclass="Chemical">natioclass="Chemical">n. For the +class="Chemical">n class="Chemical">N plots, nitrogen was added in the form of ammonium nitrate (NH4NO3) resulting in 15 g of NH4NO3 and 5.25 g of N/m2. The +P plots received 0.5 g of phosphorus in the form of 2.2 g of potassium dihydrogen phosphate (KH2PO4). +NP plots received 15 g of NH4NO3 and 2.2 g of KH2PO4. For controls, stream water from the same source was sprayed onto the plots. These levels of nutrient addition were designed to result in a pulse of nutrients that would greatly overcome any possible natural limitations.

Plots were sampled prior to the application of fertilization treatment. In each plot surface soil was collected (0–5 cm) from 2 locations, and samples were composited to generate one sample per plot. Samples were obtained in the same manner one year following the fertilization treatment. Ethanol and paper towels were used to sterilize the tools before sampling each individual plot.

Samples were collected in a similar manner along three transects of varying age across the glacial forefield both years; molecular analyses were done on the samples collected in 2011. These reference soils represented advancing stages of succession: soils that had been exposed for approximately 5 years, soils with biological soil crust formation (approximately 20 years after exposure), and soils with 25–50% vegetation cover (approximately 85 years after exposure). At the field site, samples were kept in a cooler on ice for transport to Boulder, CO. Soils were sieved (to 2 mm), and then stored at 4°C for soil characterization. A subsample was immediately archived in a −80°C freezer for molecular analysis and later used for n class="Chemical">KCl extractioclass="Chemical">ns.

Soil Analysis

Gravimetric soil moisture and pH (using a ratio of 2 g soil to 4 mL DI n class="Chemical">H2O) were assayed based oclass="Chemical">n staclass="Chemical">ndard methods [9]. For total orgaclass="Chemical">nic C aclass="Chemical">nalysis, class="Chemical">n class="Chemical">carbonate (inorganic C) removal was first performed on dried, ground soils [9]. 50 mg of these processed soils were packed into tin capsules; %C and %N were determined using a Thermo Finnigan EA 1112 Series Flash Elemental Analyzer (Thermo Fisher Scientific, Inc., Waltham, Massachusetts, USA) [31]. Bio-available P concentrations were measured on air-dried and sieved soil (2 mm×2 mm) by extracting 3–5 g of soil with 0.5 M sodium bicarbonate for 30 minutes [32]. Extracts were filtered and analyzed colorimetrically using the ammonium molybdate-malachite green method [33] adapted for microplate analysis. NH4 + and NO3 −/NO2 − extractable N were analyzed from soils using 2M KCl with 1 hour shaking and a 22 hour extraction period [34]. This analysis was performed on soils that were frozen at −80°C. Although not fresh samples, these soils typically withstand extreme fluctuations in temperature [35] and the data presented here are intended for within study comparison only. NH4 + and NO3 −/NO2 − were measured on a Lachat QuikChem 8500 Flow Injection Analyzer (Lachat Instruments, Hach Company, Loveland, CO) and BioTek Synergy 2 Multidetection Microplate Reader (BioTek, Winooski, VT) respectively.

DNA Extractions for 454 pyrosequencing

MO BIO PowerSoil™ Dn class="Chemical">NA Isolatioclass="Chemical">n kits were used as per the maclass="Chemical">nufacturer's iclass="Chemical">nstructioclass="Chemical">ns for Dclass="Chemical">n class="Chemical">NA extractions of soil samples (Mo Bio Laboratories, Inc., Carlsbad, CA). PCR-amplified bacterial 16S rRNA genes from the genomic DNA of the soil samples were generated using a universal bacterial 27F and 338R primer set as described by Hamady et al. [36], and reaction conditions followed those described by Fierer et al. [37], though modified to 25 PCR cycles. Primers included a 2 bp linker, the 454 Roche Titanium A/B primer, and a unique, 12 base pair error-correcting Golay barcode for pyrosequencing as detailed by Knelman et al. [21]. 454 Life Sciences GS FLX Titanium pyrosequencing of the 16S rRNA gene amplicons was completed by the Duke Institute for Genome Sciences & Policy (Duke University, North Carolina).

Pyrosequence and statistical analysis

Using QIIME, sequences were limited to those of a sequence length of 200 to 400 base pairs, a maximum of 5 homopolymers, a minimum quality score of 25, and a maximum of ambiguous bases/primer mismatches of 0; reverse primers were removed, and all samples were then denoised using flowgram clustering in QIIME [38]. Chloroplast sequences were removed. OTUs were selected at a 97% identity level by clustering based on representative sequences via UCLUST [39]. The Ribosomal Database Classifier [40], a naïve Bayesian classifier, was employed to assign taxonomic identification to OTUs. After sequence alignments based on the n class="Chemical">NAST algorithm [41], a phylogeclass="Chemical">ny was coclass="Chemical">nstructed with the FastTree algorithm [42]. OTU tables were rarified to the lowest class="Chemical">number of sequeclass="Chemical">nces iclass="Chemical">n a sample: 407 for commuclass="Chemical">nity dissimilarity aclass="Chemical">nalyses of fertilizatioclass="Chemical">n plots. Refereclass="Chemical">nce traclass="Chemical">nsects of advaclass="Chemical">nciclass="Chemical">ng age iclass="Chemical">ncluded 6, 5, aclass="Chemical">nd 3 sequeclass="Chemical">nced replicate samples, respectively, aclass="Chemical">nd were rarefied to 71 to iclass="Chemical">nclude all of these samples. For comparisoclass="Chemical">n of refereclass="Chemical">nce samples aclass="Chemical">nd fertilizatioclass="Chemical">n plots this workflow was repeated. Iclass="Chemical">n order to examiclass="Chemical">ne differeclass="Chemical">nces amoclass="Chemical">ng bacterial commuclass="Chemical">nities, pairwise distaclass="Chemical">nce matrices based oclass="Chemical">n weighted Uclass="Chemical">niFrac, a phylogeclass="Chemical">netic distaclass="Chemical">nce metric, were geclass="Chemical">nerated for eclass="Chemical">ntire commuclass="Chemical">nities aclass="Chemical">nd the cyaclass="Chemical">nobacterial subset of commuclass="Chemical">nities iclass="Chemical">n fertilizatioclass="Chemical">n plots [43], [44]. The Priclass="Chemical">ncipal Coordiclass="Chemical">nate Aclass="Chemical">nalysis (PCoA) ordiclass="Chemical">natioclass="Chemical">ns were coclass="Chemical">nstructed based oclass="Chemical">n OTU tables aclass="Chemical">nd weighted Uclass="Chemical">niFrac distaclass="Chemical">nce matrices for overall commuclass="Chemical">nities. The QIIME-geclass="Chemical">nerated OTU tables were used to evaluate the relative abuclass="Chemical">ndaclass="Chemical">nce of all taxa.

Primer v6 software [45] was used to perform permutational An class="Chemical">NOVAs (PERMAclass="Chemical">n class="Chemical">NOVA) to compare phylogenetic distances among bacterial communities. PERMANOVA tests were used on both UniFrac beta diversity matrices of the entire communities and cyanobacterial portions of communities. PERMANOVA analysis was also employed to assess differences among treatment-affected communities and successional reference communities. For all comparisons with reference communities, data were rarefied to the lowest sampling depth among both fertilization plot and reference plot samples.

R software [46] was used for further statistical analysis. The PERMDISP2 procedure (with permutational P-values) from the R vegan package to test homogeneity of group dispersions (variances) was also employed via QIIME in order to test for differences in community phylogenetic dispersion (UniFrac) in fertilized samples and reference successional communities [47], [48]. As well, the pgirmess package in R was used to evaluate comparisons among reference chronosequence soil relative abundance data via the Kruskal Wallis test. To assess treatment vs. temporal effects underlying shifts in overall phylogenetic community composition, a Tukey's n class="Disease">HSD post-hoc test was used to compare Uclass="Chemical">niFrac distaclass="Chemical">nces of paired pre- aclass="Chemical">nd post-treatmeclass="Chemical">nt plots with paired coclass="Chemical">ntrol plots from both years. Additioclass="Chemical">nally, to assess the relative abuclass="Chemical">ndaclass="Chemical">nces of bacterial taxa, we compared the differeclass="Chemical">nces iclass="Chemical">n paired pre- to post-treatmeclass="Chemical">nt taxoclass="Chemical">n relative abuclass="Chemical">ndaclass="Chemical">nces for each treatmeclass="Chemical">nt with that of paired coclass="Chemical">ntrol plots via Tukey's class="Chemical">n class="Disease">HSD post-hoc tests. To examine the relationship between treatment-related community shifts from our fertilization experiment and reference communities across advancing stages of soil development, we examined the relationship between weighted UniFrac phylogenetic dissimilarity and time between +NP communities and reference communities via a Spearman correlation Mantel test. The Mantel test tests the null hypothesis that there is no correlation between +NP and reference community dissimilarity and chronosequence age rank.

All relative abundance data and environmental variables were evaluated for normality. Taxon relative abundances and fertilizer plot n class="Chemical">NO3 −/class="Chemical">n class="Chemical">NO2 − were square root transformed to achieve a normal distribution prior to statistical analysis. All other edaphic factors were natural log transformed. ANOVAs, Tukey's HSD, and Kruskal Wallis post-hoc tests were used to assess differences in pH, %C, P, N pools and soil moisture in fertilization plots and reference chronosequence soils. Percent N was below the detection limit in a majority of samples and thus removed from statistical evaluations.

Sequences and metadata have been deposited in Ficlass="Disease">gShare aclass="Chemical">nd are available with the DOIs: 10.6084/m9.ficlass="Chemical">n class="Disease">gshare.1050042 (metadata) and 10.6084/m9.figshare.1048992 (sequences).

Results and Discussion

Together, our analyses demonstrate that a single +n class="Chemical">NP applicatioclass="Chemical">n caused the bacterial commuclass="Chemical">nity structure of the 3-year-old barreclass="Chemical">n soils to coclass="Chemical">nverge with the structure of 85-year-old vegetated soils after oclass="Chemical">nly oclass="Chemical">ne year. First, paired pre- aclass="Chemical">nd post-treatmeclass="Chemical">nt plot commuclass="Chemical">nity differeclass="Chemical">nces (weighted Uclass="Chemical">niFrac distaclass="Chemical">nce) were assessed amoclass="Chemical">ng all plot categories usiclass="Chemical">ng aclass="Chemical">n Aclass="Chemical">n class="Chemical">NOVA. The +NP plots showed a significant community shift in response to the treatment (Tukey's HSD; P = 0.037); no other significant differences in community structure were detected between treatments and controls (Tukey's HSD; P>0.05). A PCoA ordination (Fig. 1) revealed a successional trend in community composition across the reference chronosequence, with post-treatment +NP communities clustering with the oldest reference communities. A PERMANOVA analysis demonstrated that there were no significant differences among pre-treatment communities (Table 1). However, communities in post-treatment +NP plots were significantly different from both pre- and post-treatment controls, including the paired pre-treatment +NP plots (PERMANOVA, P<0.05, Table 1). When +NP communities were compared to reference communities across the natural chronosequence, a Mantel test of pairwise average UniFrac [43], [44] distances between +NP plots and reference samples revealed significant patterns of decreasing dissimilarity: +NP communities were most similar to the 85 year old successional soils (Fig. 2, ρM = −0.35 P = 0.01). The PERMANOVA analysis also showed that +NP communities were significantly different than communities of all successional stages except those of the oldest transect (85 years old) (Table 1). These results suggest that fertilization drives community composition away from early successional stages and results in convergence with communities of older soils. Likewise, the phylogenetic dispersion [47], [48] of +NP communities was significantly different from all reference communities except those in the 85 year old soils (Table 2). We note that our PERMANOVA analysis was not corrected for multiple comparisons due to the low statistical power of our study, but the general results of this analysis were nonetheless corroborated by our other statistical analyses of treatment effect (ANOVA/Tukey's HSD of pre- and post-treatment community shifts) and convergence of the +NP plots to the oldest successional soils (Mantel test of +NP community distance compared to successional reference samples over time).

Figure 1

Principal Coordinates Analysis (PCoA) ordination plot of bacterial communities from the field fertilization experiment and bacterial communities from the successional chronosequence.

Only the +NP treatment communities are shown because the +N and +P treatments did not result in significant community shifts. PCoA visually represents differences among community composition as the distance between points. Triangles represent communities from the natural chronosequence: red  = 5 years old; orange  = 20 years old; blue  = 85 years old. Circles represent communities from the fertilization experiment: black  =  pre-treatment control; grey  =  post-treatment control; purple  =  pretreatment +NP; Pink  =  post-treatment +NP. Our analysis revealed significant community shifts over the reference chronosequence (triangles) as well as a significant response to +NP fertilization (circles). As well, the PCoA analysis demonstrates that the +NP communities (pink circles) group with the oldest soils from the chronosequence (blue triangles).

Table 1

Post-treatment +NP phylogenetic community structure was significantly different from controls and from all communities from the reference chronosequence with the exception of communities in the oldest soils (P<0.05).

Permutational MANOVA (PERMANOVA) contrast P-values
Sample vs. Sample	post-treatment control	post-treatment +N	post-treatment +P	post-treatment +NP
pre-treatment control	0.415	0.066	0.031	0.026
post-treatment control	---	0.422	0.072	0.036
pre-treatment +NP-paired	0.114	0.085	0.024	0.032
post-treatment +NP plots	0.036	0.023	0.18	---
succession timepoint 1	0.124	0.003	0.005	0.006
succession timepoint 2	0.105	0.009	0.018	0.022
succession timepoint 3	0.152	0.055	0.179	0.162
Significant P-values (P<0.05) bolded.

Controls showed no differences from any contrasts (P>0.05). Significant P-values (P<0.05) are bolded.

Figure 2

Relationship between +NP treatment-affected communities and reference communities.

A box plot shows the average weighted UniFrac [43], [44] distance between +NP-treated communities and reference communities with increasing successional time. A Mantel test demonstrates that +NP communities show decreasing dissimiliarty as compared to the reference communities over advancing stages of succession (ρM = −0.35 P = 0.01).

Table 2

Post-treatment +NP communities showed differences from all reference succession communities with the exception of the oldest transect (P<0.05).

Homogeneity of Dispersion (PERMDISP) P-values
Sample vs. Sample	post-treatment +NP	post-treatment control
succession timepoint 1	0.022	0.508
succession timepoint 2	0.042	0.588
succession timepoint 3	0.555	0.997
Significant P-values (P<0.05) bolded.

Controls showed no difference in community dispersion from communities of any of the reference succession transects (P>0.05). Significant P-values (P<0.05) are bolded.

Only the +n class="Chemical">NP treatmeclass="Chemical">nt commuclass="Chemical">nities are showclass="Chemical">n because the +class="Chemical">n class="Chemical">N and +P treatments did not result in significant community shifts. PCoA visually represents differences among community composition as the distance between points. Triangles represent communities from the natural chronosequence: red  = 5 years old; orange  = 20 years old; blue  = 85 years old. Circles represent communities from the fertilization experiment: black  =  pre-treatment control; grey  =  post-treatment control; purple  =  pretreatment +NP; Pink  =  post-treatment +NP. Our analysis revealed significant community shifts over the reference chronosequence (triangles) as well as a significant response to +NP fertilization (circles). As well, the PCoA analysis demonstrates that the +NP communities (pink circles) group with the oldest soils from the chronosequence (blue triangles).

A box plot shows the average weighted UniFrac [43], [44] distance between +n class="Chemical">NP-treated commuclass="Chemical">nities aclass="Chemical">nd refereclass="Chemical">nce commuclass="Chemical">nities with iclass="Chemical">ncreasiclass="Chemical">ng successioclass="Chemical">nal time. A Maclass="Chemical">ntel test democlass="Chemical">nstrates that +class="Chemical">n class="Chemical">NP communities show decreasing dissimiliarty as compared to the reference communities over advancing stages of succession (ρM = −0.35 P = 0.01).

Our results suggest that nutrient colimitation is an important control on microbial primary succession in this system. Because of low statistical power, it is difficult to discern whether this colimitation is simultaneous, meaning that both nutrients need to be present for a community response, or independent, meaning that each nutrient in isolation may elicit some response [49]. However, there is some evidence that single nutrient additions may cause a smaller response than when both nutrients are abundant. For example, our results show that post-treatment +P communities are not significantly different from post treatment +n class="Chemical">NP commuclass="Chemical">nities (Table 1). As well, both +class="Chemical">n class="Chemical">N and +P plots show patterns of convergence similar to +NP plots in comparison with ongoing natural succession; by contrast, control plots do not display convergence (Table 1). Thus, +N and +P communities may represent intermediate states between control and +NP plots, but we were not able to statistically demonstrate an underlying treatment effect.

While our study is unique as we established and resampled nutrient addition plots in a remote glacial forefield, the rapidly changing nature of the Puca Glacier landscape and criteria for setn class="Chemical">ting up plots oclass="Chemical">n a stable aclass="Chemical">nd relatively homogeclass="Chemical">nous surface limited replicatioclass="Chemical">n aclass="Chemical">nd class="Chemical">necessitated rarefactioclass="Chemical">n of sequeclass="Chemical">nciclass="Chemical">ng depth to 71 to iclass="Chemical">nclude all available samples. As such, we ackclass="Chemical">nowledge the class="Chemical">need to be circumspect iclass="Chemical">n drawiclass="Chemical">ng coclass="Chemical">nclusioclass="Chemical">ns as such factors curbed the statistical power of our study aclass="Chemical">nd poteclass="Chemical">ntially our ability to detect smaller magclass="Chemical">nitude treatmeclass="Chemical">nt effects iclass="Chemical">n the +class="Chemical">n class="Chemical">N and +P additions, for example. However, we note that the patterns shown here are robust to even lower rarefaction depths (55–70); thus, it is likely that observed patterns are real. Nonetheless, our research shows the greatest, and only statistically significant treatment effect on microbial communities under +NP additions, suggesting the effect of both nutrients in tandem is important in succession.

Interesn class="Chemical">tingly, staclass="Chemical">ndiclass="Chemical">ng class="Chemical">nutrieclass="Chemical">nt pool aclass="Chemical">nalysis leclass="Chemical">nds some iclass="Chemical">nsight iclass="Chemical">nto particular dyclass="Chemical">namics that may uclass="Chemical">nderlie class="Chemical">nutrieclass="Chemical">nt colimitatioclass="Chemical">n iclass="Chemical">n this autotrophic chroclass="Chemical">nosequeclass="Chemical">nce. For example, +P aclass="Chemical">nd +class="Chemical">n class="Chemical">NP soils both show significant increases in ammonium pools in comparison with control plot soils (Table S1), which is consistent with a body of research that demonstrates P limitation is a strong control of N-fixation [50], [51], and may be particularly strong in this autotrophic chronosequence that features cyanobacterial N-fixers [9]. Likewise, +N plots show a significant increase in bioavailable-P relative to control plots (Table S1), a pattern supported by research that shows N availability may limit the production of phosphatase enzymes [51]–[53]. Thus, these particular biochemical pathways lead to a coupling of nutrient cycles, which appears to be reflected in a colimitation to successional processes.

Despite the multitude of well documented changes across successional gradients including shifts in pH, C pools, plant cover and biotic historical factors, nutrient addition alone not only caused changes in early successional community structure, but induced convergence with late successional soil communities of the natural chronosequence (Fig. 1 and 2 and Tables 1 and 2). For example, strong changes in %C, another known filter on microbial communities, were observed across the natural chronosequence but not in +n class="Chemical">NP plots (Tables S1 aclass="Chemical">nd S2). Iclass="Chemical">n other ecosystems, the effects of fertilizatioclass="Chemical">n oclass="Chemical">n microbial commuclass="Chemical">nity structure have beeclass="Chemical">n attributed to chaclass="Chemical">nges iclass="Chemical">n placlass="Chemical">nt productivity or commuclass="Chemical">nity structure [14]. However, it is importaclass="Chemical">nt to class="Chemical">note that while the +class="Chemical">n class="Chemical">NP fertilization caused sparse vegetation (<15 cm tall) to colonize after one year at our site, soils were collected at least 75 cm from these small plants. Altogether, our results suggest that the effects of the +NP fertilization on microbial community succession were direct and not mediated through changes in other aspects of the abiotic environment or through the effects of plants on soil communities.

Our field-based fertilization experiment helps to extend exisn class="Chemical">ting ecological theory regardiclass="Chemical">ng the role of class="Chemical">nutrieclass="Chemical">nt limitatioclass="Chemical">ns iclass="Chemical">n successioclass="Chemical">n [4], [5], [54] to microbial commuclass="Chemical">nities preseclass="Chemical">nt iclass="Chemical">n the earliest primary successioclass="Chemical">nal soils, which are importaclass="Chemical">nt for biogeochemical cycliclass="Chemical">ng, physical soil developmeclass="Chemical">nt, aclass="Chemical">nd placlass="Chemical">nt coloclass="Chemical">nizatioclass="Chemical">n [9], [21], [30]. While it is widely ackclass="Chemical">nowledged that microbes caclass="Chemical">n alter soil fertility aclass="Chemical">nd class="Chemical">nutrieclass="Chemical">nt cycliclass="Chemical">ng processes, aclass="Chemical">nd that chaclass="Chemical">nges iclass="Chemical">n soil class="Chemical">nutrieclass="Chemical">nt pools aclass="Chemical">nd microbial commuclass="Chemical">nities occur over primary successioclass="Chemical">n [9], [12], [13], [30], to what exteclass="Chemical">nt class="Chemical">nutrieclass="Chemical">nts directly structure soil microbial commuclass="Chemical">nities is class="Chemical">not clear. Our fertilizatioclass="Chemical">n experimeclass="Chemical">nt allowed us to decouple the effects of chaclass="Chemical">nges iclass="Chemical">n microbial commuclass="Chemical">nities oclass="Chemical">n class="Chemical">nutrieclass="Chemical">nt cycles aclass="Chemical">nd to directly democlass="Chemical">nstrate the iclass="Chemical">nflueclass="Chemical">nce of class="Chemical">nutrieclass="Chemical">nt pools oclass="Chemical">n microbial successioclass="Chemical">n. Correlative studies are less powerful because they caclass="Chemical">nclass="Chemical">not isolate the impact of iclass="Chemical">ndividual factors amidst the multiplicity of soil properties that chaclass="Chemical">nge with successioclass="Chemical">n, aclass="Chemical">nd because measured soil properties may be decoupled from microbial commuclass="Chemical">nity compositioclass="Chemical">n iclass="Chemical">n time.

Despite the high fertilization rate we used, the nutrient addition treatment did not push communities to an alternative or novel state, but simply accelerated succession, rapidly producing a community that was structurally most similar to the community in the 85 year old soils in the chronosequence (Fig. 1 and 2 and Tables 1 and 2). Thus, our data highlight the stability of soil microbial communities [55]. Few studies have explicitly evaluated nutrients in the context of longer-term successional reference plant communities to understand how nutrients may either drive succession or shape alternative stable states in communities. However, in a study of salt marsh vegetation, Van Wijnen and Bakker [56] observed that fertilization of young marsh communities resulted in plant communities that resembled those of older, unfertilized marshes. These results further suggest that nutrient-related mechanisms for succession may be generalizable between plant and microbial communities.

The relative abundance of cyanobacteria significantly increased in the +n class="Chemical">NP plots aclass="Chemical">nd the phylogeclass="Chemical">netic structure of the cyaclass="Chemical">nobacterial commuclass="Chemical">nities iclass="Chemical">n post-treatmeclass="Chemical">nt +class="Chemical">n class="Chemical">NP plots was significantly different from the paired pre-treatment +NP and pre-/post-treatment control plots (PERMANOVA, P<0.05, Table S3). Although not significant, cyanobacterial relative abundance nearly doubled between the oldest and youngest stages of the reference chronosequence and past work at this site has documented similar successional changes in cyanobacterial community structure (Table S2) [9], [30]. Consistent with these results, a laboratory experiment evaluating microbial autotrophs from this site demonstrated that P additions resulted in significant increases in the growth rate of photoautotrophic crusts [28]. Both N fixation rates and the relative abundance of N-fixing cyanobacteria show successional trends at this site as well [30], suggesting that N availability may also limit microbial growth and activity. The current study adds to this work and demonstrates that both N and P together are important colimiting controls over community successional processes in this system (Tables 1 and 2).

The increase in the relative abundance of cyanobacteria in the +n class="Chemical">NP plots may reflect their ecological advaclass="Chemical">ntage iclass="Chemical">n this low C eclass="Chemical">nviroclass="Chemical">nmeclass="Chemical">nt. Iclass="Chemical">n a laboratory study, Drakare [57] observed that P additioclass="Chemical">ns eclass="Chemical">nhaclass="Chemical">nced cyaclass="Chemical">nobacterial populatioclass="Chemical">ns, but oclass="Chemical">nly iclass="Chemical">n aclass="Chemical">n eclass="Chemical">nviroclass="Chemical">nmeclass="Chemical">nt where low C coclass="Chemical">nceclass="Chemical">ntratioclass="Chemical">ns coclass="Chemical">nstraiclass="Chemical">ned heterotrophic growth. Iclass="Chemical">ncubatioclass="Chemical">n studies of early successioclass="Chemical">nal soils that fouclass="Chemical">nd iclass="Chemical">ncreases iclass="Chemical">n heterotrophic activity iclass="Chemical">n respoclass="Chemical">nse to both class="Chemical">n class="Chemical">N and C (but not to N alone) are also consistent with this interpretation [58], [59]. These results indicate that C often limits the response of the heterotrophic community to nutrient additions, whereas cyanobacteria can readily take advantage of such nutrients to fuel photosynthesis. By extension, we argue that the observed effects of N and P additions on microbial community succession are likely to apply only to autotrophic successional sequences, and that heterotrophic succession (sensu Fierer et. al [27]) may be controlled by a different suite of resources, including C availability.

Microbes are fundamental to soil physical and chemical development and underlie ecosystem function, thus understanding the factors that drive soil microbial community succession is key to predicn class="Chemical">ting aclass="Chemical">nd maclass="Chemical">nagiclass="Chemical">ng ecosystem developmeclass="Chemical">nt. Particularly iclass="Chemical">n low class="Chemical">nutrieclass="Chemical">nt eclass="Chemical">nviroclass="Chemical">nmeclass="Chemical">nts, microbial activity has major effects oclass="Chemical">n soil, placlass="Chemical">nt commuclass="Chemical">nity, aclass="Chemical">nd ecosystem developmeclass="Chemical">nt [9], [30], [60], [61]. Likewise, low class="Chemical">nutrieclass="Chemical">nt eclass="Chemical">nviroclass="Chemical">nmeclass="Chemical">nts may feature more promiclass="Chemical">neclass="Chemical">nt class="Chemical">nutrieclass="Chemical">nt colimitatioclass="Chemical">ns [49]. As such, the results of this study have importaclass="Chemical">nt implicatioclass="Chemical">ns for uclass="Chemical">nderstaclass="Chemical">ndiclass="Chemical">ng class="Chemical">nutrieclass="Chemical">nt coclass="Chemical">ntrols oclass="Chemical">n ecosystem developmeclass="Chemical">nt aclass="Chemical">nd relevaclass="Chemical">nt models for microbial successioclass="Chemical">n. Furthermore, while early successioclass="Chemical">nal microbial commuclass="Chemical">nities may vary stroclass="Chemical">ngly iclass="Chemical">n both compositioclass="Chemical">n aclass="Chemical">nd iclass="Chemical">n terms of the specifics of resource availability (e.g., heterotrophic vs. autotrophic), our study provides evideclass="Chemical">nce that class="Chemical">nutrieclass="Chemical">nt colimitatioclass="Chemical">n may provide a geclass="Chemical">neralizable mechaclass="Chemical">nism for microbial commuclass="Chemical">nity successioclass="Chemical">n iclass="Chemical">n autotrophic successioclass="Chemical">nal sequeclass="Chemical">nces. Our data also support receclass="Chemical">nt evideclass="Chemical">nce for the stability of soil microbial commuclass="Chemical">nities, as fertilizatioclass="Chemical">n simply accelerated successioclass="Chemical">n aclass="Chemical">nd did class="Chemical">not push commuclass="Chemical">nities iclass="Chemical">nto a class="Chemical">novel state. Overall, the details of microbial class="Chemical">nutrieclass="Chemical">nt limitatioclass="Chemical">ns preseclass="Chemical">nted hereiclass="Chemical">n are esseclass="Chemical">ntial to uclass="Chemical">nderstaclass="Chemical">ndiclass="Chemical">ng the factors that structure early successioclass="Chemical">nal microbial commuclass="Chemical">nities, the profouclass="Chemical">nd coclass="Chemical">ntributioclass="Chemical">ns they make to soil developmeclass="Chemical">nt, aclass="Chemical">nd the ecosystem processes they mediate.

Mean of Edaphic Properties and Tukey's n class="Disease">HSD Comparisoclass="Chemical">ns for Post-Treatmeclass="Chemical">nt Plots.

(DOCX)

Click here for additional data file.

Mean of Edaphic Properties and Cyanobacterial Relative Abundance for Reference Chronosequence.

(DOCX)

Click here for additional data file.

Mean Relative Abundances of Major Taxa at Puca Glacier Site.

(DOCX)

Click here for additional data file.