Phosphorus, not nitrogen, limits plants and microbial primary producers following glacial retreat

See allHide authors and affiliations

Science Advances  23 May 2018:
Vol. 4, no. 5, eaaq0942
DOI: 10.1126/sciadv.aaq0942


Current models of ecosystem development hold that low nitrogen availability limits the earliest stages of primary succession, but these models were developed from studies conducted in areas with temperate or wet climates. Global warming is now causing rapid glacial retreat even in inland areas with cold, dry climates, areas where ecological succession has not been adequately studied. We combine field and microcosm studies of both plant and microbial primary producers and show that phosphorus, not nitrogen, is the nutrient most limiting to the earliest stages of primary succession along glacial chronosequences in the Central Andes and central Alaska. We also show that phosphorus addition greatly accelerates the rate of succession for plants and for microbial phototrophs, even at the most extreme deglaciating site at over 5000 meters above sea level in the Andes of arid southern Peru. These results challenge the idea that nitrogen availability and a severe climate limit the rate of plant and microbial succession in cold-arid regions and will inform conservation efforts to mitigate the effects of global change on these fragile and threatened ecosystems.


A long-held paradigm of biogeochemical and ecological theory is that early successional ecosystems are primarily nitrogen (N)–limited, whereas later successional ecosystems are primarily phosphorus (P)–limited (14). This hypothesis grew out of classic ecological studies and the observation that all the elements required for plant growth, except N, are found in geologic substrates on which primary succession occurs and that soils tend to weather from young, P-rich systems to older, P-poor systems (1, 5). However, this paradigm was developed and has mostly been tested in humid coastal and island ecosystems (1, 58) and, to a more limited extent, in semiarid temperate ecosystems (2, 9). Although some studies have shown that early successional systems can be P-limited (7), to date, this theory has not been tested in cold-arid inland environments where glacial and ice cap retreat due to global warming is opening up vast new landscapes to primary succession. These environments are among the systems most at risk because of global warming (10) and are also areas where low temperatures and aridity may further limit the rate of ecosystem succession beyond limitations imposed by low nutrient availability. It has been shown that the preplant stages of ecosystem succession are protracted in these environments (11, 12) and, in the most extreme cold-arid ecosystems, succession never progresses beyond the stage where microbes are the only primary producers of the system (13, 14). Therefore, it is important to test models of ecosystem succession for both microbes and plants in cold-arid ecosystems, especially in areas far from the moderating influences of coastal climates.

Here, we show that microbial and plant succession are primarily P-limited at cold inland sites in both the Alaska Range and High Andes of Peru. We used in situ fertilization experiments and realistic microcosms to show that phototrophic organisms ranging in size from a few micrometers in diameter to macroscopic plants are primarily P-limited. Our field nutrient addition experiments conducted at high-elevation sites [5000 m above sea level (m.a.s.l.)] show that phosphorus limitation overrides constraints of a cold-arid climate and that both microbial and macroscopic primary producers are primarily P-limited during the early stages of primary succession. These experimental results are supported by evidence from biogeochemical stoichiometry and laboratory incubations of soils from early successional soils in the Canadian and Colorado Rockies (15, 16) and from earlier studies of the Peruvian and inland Alaskan sites discussed in the present study (12, 1719).


Our primary research site, the rapidly deglaciating Sibinacocha watershed in arid southern Peru, has been the focus of a multinational effort to understand the biological and societal impacts of rapid deglaciation (11, 1721). Previous work indicated that primary succession at this site is extremely slow, with limited plant colonization even 150 years after glacial retreat. This slow rate of plant colonization has been attributed to N limitation and to the extreme climatic conditions at this site, that is, aridity for much of the year combined with extreme temperature fluctuations due to the elevation (5000 to 6000 m.a.s.l.) of the watershed (11). To test the hypothesis that N is limiting at this site, we established a long-term, full-factorial N × P fertilization study in 2010. Nutrients (5.25 g of N as NH4NO3 and/or 0.50 g of P as KH2PO4 per square meter) were added only once at the beginning of the experiment. It should be noted that the +P treatment also included K, but it is unlikely that the addition of this small amount of K would have an effect in this system because the soil is already rich in K; the total pool of K in these early successional soil is about 40 times larger than the total pool of P (table S1). Despite the one-time pulse of nutrients, plants showed a strong response to nutrient addition over the next 6 years, with an exponential increase in the number of plants per plot in treatments receiving K2HPO4 (henceforth called P), whereas plots that received just N were barely distinguishable from controls even after 6 years (Fig. 1A). When plant colonization was high enough to measure using percent cover approaches [point intercept (p-i) and modified normalized difference vegetation index (mNDVI)], these more robust analyses also revealed P limitation (Figs. 1B and 2). Plots receiving P and P+N had high percent coverage of the soil surface, whereas control and +N treatments had minimal to no coverage. Not all plant species could be identified because most plants were very small and lacked reproductive structures, but those that could be identified were Valeriana sp. (either Valeriana pycnantha or Valeriana nivalis), Arenaria sp. (likely Arenaria digyna), and Deyeuxia sp. (likely Deyeuxia hackelii; visible in Fig. 2C). Our results show a very clear signal of P limitation in these plants as opposed to colimitation by N and P (22, 23). This conclusion is also reflected in the effect sizes of nutrient additions; P explained 78% of the variation in percent cover (p-i) in 2015 and 79% in 2016 (P < 0.00001, ANOVA), but N explained less than 4% in both years (P > 0.43). In addition, mNDVI data showed almost identical trends to the p-i and stem counts (Fig. 1) but this technique was also sensitive enough to detect some microbial phototrophs even in the control plots (Fig. 2) as has also been shown using NDVI in other arid systems (24).

Fig. 1 Three different measures of plant growth in response to a one-time addition of nutrients to soils (in August 2010) that were uncovered by the Puca Glacier 3 years preceding the start of the experiment.

This site was initially barren, but plots with added P exhibited exponential increase in the number of plants within (A). By the fifth year (2015), plant biomass was high enough that it was feasible to measure the percent cover by plants within the plots. This was done manually using the p-i method and digitally using mNDVI in 2015 and 2016 (B). The p-i method yielded almost identical results for 2015 (left) and 2016 (right), and analysis of variance (ANOVA) showed an overwhelmingly large effect size of P addition (ω2, inset, pie charts). mNDVI analysis showed very similar results (bottom of Fig. 3B), although the effect size of P was smaller (***P < 0.001, **P < 0.01). The effect of N addition was not statistically significant in any of these four analyses.

Fig. 2 Near-infrared–enhanced photographs of field plots and mNDVI comparison.

In these false-color photographs, the visible light channels of the sensor have been averaged and rendered as a grayscale background, and the near-infrared (NIR) channel has been overlayed in red. Control plots (the dashed line is the plot boundary) (A) did not contain any visible plants, but some pixels in the plot area still had mNDVI values above the threshold of 0.1 (B), possibly indicating detection of microbial phototrophs. However, plots with P added (C) contained many plants, which were detected with the NIR sensor and statistically differentiated from barren areas using mNDVI. The mNDVI distribution plot (B) shows that the +P plot has many more pixels above the mNDVI threshold (dashed red line) than the control plot.

Although mNDVI robustly quantified the influence of P and N on plants and detected microbial phototrophs even in the control plots (Figs. 1B and 2), our sensor had insufficient optical resolution to quantify microbial responses to P and/or N. To overcome this problem, we used microscopic counts of phototrophic colonies and chlorophyll autofluorescence microscopy to quantify the effects of nutrients on microbial phototrophic biomass and DNA sequencing to determine how nutrients change the structure of the phototrophic community (Fig. 3). DNA sequencing revealed that the relative abundance of cyanobacteria in the order Nostocales increased significantly (P < 0.00005, ANOVA) in response to added P, with 79% of the variance in Nostocales relative abundance explained by P; in contrast, other cyanobacteria showed no significant response to P or N (Fig. 3A). The Nostocales contain many N fixers, and a metagenome of the dominant phylotype showed an abundance of nifH genes from the genus Nostoc. The increase in relative abundance of N-fixing cyanobacteria after addition of P (Fig. 3A) is consistent with P being the ultimate limiting nutrient in this system as defined by Vitousek et al. (25). In addition, relatively large (up to 300 μm in diameter) spherical Nostoc colonies growing on the soil surface (Fig. 3, B and C) were significantly more abundant in plots with added P (Fig. 3D). We also used chlorophyll autofluorescence microscopy (Fig. 3E) to quantify even smaller (down to the 1- to 10-μm range) microbial phototrophs, and these measurements also demonstrated a significant response to P (Fig. 3F). Both plants and microbes responded positively to N as well, possibly indicating slight simultaneous colimitation (22, 23, 26) or community colimitation (27) where microbial N fixers are principally P-limited, thereby alleviating N limitation or NP colimitation for plants. However, the effect size of N in our field experiment is small (Figs. 1 and 2), and it is clear that P limits the system as a whole in regard to primary productivity. Together, our results show that organisms ranging in size from a few micrometers up to plants that are five orders of magnitude larger are all primarily P-limited in this early successional ecosystem.

Fig. 3 Shifts in cyanobacterial communities from nutrient addition and nutrient limitation of microbial phototrophs.

Both N and P addition changed cyanobacterial community structure, but P addition caused a marked increase in the relative abundance of N-fixing Nostocales (pink) (A). Nostoc were visually apparent across several orders of magnitude in samples from the field plots as large spherical colonies (B) or strands of cells (C) and were identified using shotgun-metagenomic and amplicon sequencing as Nostoc commune (16S ribosomal RNA gene was 100% identical to GenBank accession no. KY380004). Microscopic counts of spherical Nostoc colonies showed P limitation (D), as did microscopic percent cover based on chlorophyll autofluorescence (E and F), similar to the results for macroscopic plants (Figs. 1 and 2) (***P < 0.001, **P < 0.01).

To further test P limitation at this site and to link these findings to sites in central Alaska, we used a laboratory bioassay approach to distinguish between P and/or N limitation of soil phototrophs (28, 29). At both sites, P addition resulted in higher percent cover by microbial phototrophs (Fig. 4A), and the effect size of P addition was much larger than that of N addition (Fig. 4B). Although N still had a statistically significant effect on primary productivity in both cases, the difference in effect size between N and P makes it clear that these microcosms were P-limited (22, 23). In addition, residual N was found only in the +N treatments of both experiments (Fig. 4C), indicating that P addition allowed for use of added N, but N added to the +N microcosms could not be used without the added P (as in the +N+P microcosms). On the basis of work in wetter systems, it might be thought that the N that remains unused because of P limitation would be lost because of leaching in the field. However, in arid systems, nitrate leaching is less of an issue, and nitrate levels were almost an order of magnitude higher in the plots receiving N compared to controls after 1 year in the field experiment reported here (17).

Together, and with support from other studies (11, 12, 15, 1719), these results point to a new paradigm for primary succession in cold inland ecosystems compared to well-studied coastal and island landscapes. Despite high concentrations of P in the primary minerals on which succession occurs (18), P availability is highly restricted in climates where rock weathering is inhibited by cold temperatures and aridity. These conditions result in the severe P limitation described and explain the slow pace of succession of plant and microbial primary producers in places like the High Andes, Central Alaska, and perhaps in other understudied inland ranges such as the Himalaya and Karakoram.

Fig. 4 Results of laboratory bioassays.

Laboratory bioassay experiments (28, 29) showed P limitation of early successional soils from the High Andes (left) and from the receding Toklat Glacier in central Alaska (right). The percent cover of microbial phototrophs over time is significantly higher in treatments with P added and is lower when P is not added (A). Repeated-measures ANOVA shows a much larger effect size of P than of N as well (B) (***P < 0.001, **P < 0.01). Furthermore, when N is added in the absence of added P, that N was not assimilated, indicating that the microcosms were P-limited even in their ability to immobilize N (C). FOV, fields of view [see Darcy and Schmidt (29)].


Field plots were created in August 2010 in soil at the terminus of the retreating Puca Glacier, in the Cordillera Vilcanota of Peru, at 5000 m.a.s.l. The site used was exposed in the previous year by the glacier’s retreat and faces east. Between 2010 when the plots were created and 2016 when they were last examined, the Puca Glacier retreated at an average of 28 m per year. Moraines at this site were composed of rock with high quartz and calcite content, and our experimental sites were composed mainly of glacial till with a high content of carbonates (11), although larger particles and rocks were present in the plots as well (Fig. 2). Nutrient deposition has not been quantified at this site, but at the nearby Quelccaya ice cap (30 km SE), N deposition has been thought to derive mainly from Amazonian sources in the wet season and from lightning in the dry season (30); pollen concentrations on the Quelccaya ice cap are among the highest ever observed (31). Mean annual temperature at the Puca Glacier site is 0°C and mean annual precipitation is 702 mm (average over 15 years) (21), but much of this precipitation is snow, which sublimates back to the atmosphere; thus, the soils are usually very dry when not covered in snow (11, 1618). A map of the experimental plots can be seen in fig. S1. These 16 permanent plots (1 × 1 m2) were marked using nails and a string, with small aluminum numbered tags for future identification. Plot treatment identity was assigned randomly from the four treatments of our factorial N and P design such that each treatment identity (Control, +N, +P, and +N+P) was represented by four plots.

All plots were completely devoid of visible plants or lichens at the time of nutrient addition. Nutrients were added using spray bottles as described elsewhere (17) such that applications consisted of 5.25 g of N m−2 as NH4NO3 for +N plots and 0.5 g of P as KH2PO4 for +P plots, whereas +N+P plots received both additions combined. This treatment was applied only once, at the beginning of the experiment in 2010, because the site is cold and oligotrophic, and utilization of any added nutrients was expected to be very slow. Knelman et al. (17) showed that N persisted at the site for at least 1 year after the initial treatment. Nutrients were delivered as aqueous solutions, and water without added nutrients was sprayed on the control plots to control for potential effects of water addition. The source for the water was a glacial stream at the glacier’s terminus. Soil samples were collected before treatment was added to plots and kept frozen at −20°C. The number of plant stems (vascular plants) was counted in each plot, such that each stem represented a single plant, starting in 2010 and repeated in 2012, 2013, 2015, and 2016. Stem count data were analyzed using repeated-measures ANOVA to test the influence of N and P addition on the number of plants present in the field plots.

Soils collected from the top 2 cm of plots in 2010 were used in the nutrient addition microcosm experiment by Schmidt et al. (28), which used the same N/P factorial design as the in situ field plots. Microcosm data from Schmidt et al. (28) and Darcy and Schmidt (29) were analyzed by first reducing the Toklat Glacier (Alaska, USA) data set to only include those time points that were measured in the same way as Schmidt et al. (28) and then converting field-of-view data to percent cover using the conversion equation from Darcy and Schmidt (29). Repeated-measures ANOVA was used to test the influence of N and P addition on percent cover. Briefly, the Toklat Glacier site is located in Denali National Park and Preserve, between Fairbanks and Anchorage Alaska, USA. Soils at this site are calcareous shale soils, the mean annual temperature is 2.86°C, and the mean annual precipitation (mostly as snow) is 889 mm (temperature and precipitation data from 2009) (32). Atmospheric nitrogen deposition in Denali National Park and Preserve is relatively low (33). More description of the Toklat Glacier site can be found elsewhere (29, 3436).

Soil samples were collected from each of the Puca Glacier field plots in November 2012 to analyze the composition of developing microbial communities in the soil. Soil samples were kept at −20°C until genomic DNA was extracted using a MoBio PowerSoil kit. Each extraction was amplified in triplicate using the bacterial primer set 515F/806R (37), using a unique barcode for each of the 20 samples. Each set of triplicate reactions was pooled, and amplified DNA concentrations were assayed using PicoGreen fluorimetry on a BioTek Synergy 2 microplate reader (BioTek). DNA was then diluted to equimolar concentration, pooled, and sequenced using an Illumina MiSeq (Illumina Inc.) with 2 × 150–base pair (bp) paired-end reads. 16S ribosomal DNA data were demultiplexed and processed using QIIME software (version 1.9) (38) using default settings for operational taxonomic unit (OTU) clustering and taxonomy assignment. Sequence data were rarefied to 20,000 sequences per sample. Genomes of cyanobacteria most closely related to our most abundant cyanobacterial OTU were found in the National Center for Biotechnology Information (NCBI) genome database using nucleotide BLAST (39).

Shotgun metagenomic DNA sequencing was used to determine whether the dominant cyanobacteria within +P plots had the genes necessary for N fixation. DNA from a +P plot extraction (above) was sent to Pennsylvania State University’s Genomics Core Facility, where the extracted DNA was fragmented to 400-bp lengths. Fragments were prepared for sequencing using a KAPA HyperPlus kit (KAPA Biosystems) and then sequenced using an Illumina NextSeq (Illumina Inc.) with 2 × 150–bp paired-end reads. Raw sequence data were quality-filtered to exclude sequence pairs and trim ends with an average quality score below 37 using Sickle version 1.33 ( Filtered metagenomic reads were assembled using SPAdes version 3.10.1 (40) using the “--careful” protocol and with recommended kmer values of 12, 22, 55, and 77. To cluster contigs that originated from the same species, assembled contigs were binned according to coverage level and tetramer frequencies using MaxBin version 2.2.4 (41). BLAST (39) was used to search each bin for nitrogen fixation genes and to find the bin’s closest relative in the NCBI Genome database.

We used microscopy to estimate the density of microbial phototrophs within soil samples. The abundance of spherical Nostoc colonies within soil samples was counted in 10-g samples of soil from the field plots. Each sample was placed into a 5-cm petri dish, and then, using a dissecting microscope at ×45 magnification, we counted the number of visible spherical Nostoc colonies within each dish. ANOVA was used to test the influence of N and P addition on Nostoc colony abundance. In addition, we used chlorophyll autofluorescence microscopy to compare the density of microbial phototrophs between soils from the field plots. One gram of frozen soil from each sample was placed in a 2-ml microcentrifuge tube, and 1 ml of sterile deionized water was added to it. These tubes were incubated for 24 hours with approximately 13 hours of natural sunlight to allow microbial phototrophs to reactivate. Tubes were vortexed briefly, and then, 10 μl of liquid was pipetted from the tube onto a microscope slide. A cover slip was added, and then, a chlorophyll autofluorescence lamp was used to visualize microbial phototrophs containing chlorophyll. Images were captured using a Nikon D3100 DSLR camera (Nikon) fitted to an Axioskop 2 (Carl Zeiss AG) light microscope using a 3D-printed Nikon F-mount adapter (available for download at Ten photographs of each slide were taken from different fields of view, at ×200 magnification. These photographs were analyzed in R 3.3.3 using custom software. Briefly, each image has red, green, and blue channels. Because chlorophyll autofluoresces red, only the red channel was used. When the digital photograph is read into R, the red channel is represented as a matrix of numeric values, with the dimensions of the camera’s image resolution (4000 × 3000 pixels). To determine whether a given pixel represented detection of a microbial phototroph, that pixel’s red channel intensity was compared to a threshold value. That threshold was determined by varying the threshold of a positive-control image until known phototrophs were distinguished. Pixels with red intensities over the threshold were counted as “hits,” and the percent cover (on the slide) of microbial phototrophs was calculated as hits/total × 100. ANOVA was used to test the influence of N and P addition on percent cover of microbial phototrophs.

Total mineral phosphorus (P) and potassium (K) for soil samples near the Puca Glacier’s terminus (n = 3) were determined in 2009 (18) using element analyses with a Phillips PW1400 Wavelength Dispersive Spectrometer, an x-ray fluorescence instrument. Running conditions were 40 kV and 20 mA. Samples were dried and ground to <70 μm and then bound using corn starch and pressed into a pellet. Quantitative analyses using U.S. Geological Survey rock standard GSP-1 was performed using the fundamental parameter correction function of Phillips X40 version 4.0A software.

We revisited field plots in 2015 and 2016 to take precise measurements of percent cover by plants. These measurements were taken during the wet season (February to April) so that plant cover could be recorded when it was still green. Percent cover was measured using the p-i method, which was done by laying marked strings over the plots and counting how many times the marks intersected plants. One hundred intersections were surveyed for each of the 16 plots, for both years. However, this method may miss small plants and microbial phototrophs; hence, we also used a specially made camera that can capture NIR light to perform mNDVI analysis (42). The camera was a Hero4+ (GoPro), modified to capture NIR light by the company Agribotix LLC. This modification replaced the camera’s stock filter (which blocks NIR light) with a filter that blocks red light and allows NIR light to pass through and register as red light on the camera’s sensor. Photographs of each site were taken, making sure to capture all four corners of the site using the camera’s rear screen preview. Adobe Photoshop CC version 2016 (Adobe Inc.) was used to make masks for each photograph, which are images where the plot space is colored white (R, 255; G, 255; B, 255) and everything else captured in the photograph is black (R, 0; G, 0; B, 0). The original photographs and their masks were processed in R using custom software. The mNDVI equation (an equation identical to NDVI except that the NIR channel is normalized using only one visible light channel) was used to analyze the original photograph. Previous mNDVI implementations have used red as the visible channel (42), but because our modified camera uses the red channel for NIR light, we used the blue channel instead to normalize NIR per the manufacturer’s instructions. Orbital imaging platforms (on satellites) have separate NIR and visible light sensors, which is a viable solution for multispectral imaging because the two sensors are so far from the target that the two sensors produce an image of the same focal area. However, on the ground, less than a meter from the target, the two-sensor solution is not viable (mathematically or economically). Thus, we used mNDVI instead and applied that function using the entire original photograph’s NIR channel and blue channel. Pixels that were not within the plot area were discarded, using the mask. This process was done for each plot across both years. Percent cover estimates using mNDVI and using the p-i method were tested with ANOVA.


Supplementary material for this article is available at

table S1. Element analysis results from x-ray fluorescence analysis.

fig. S1. Aerial view of the Puca Glacier field site, showing plot layout and orientation.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank R. I. Meneses, K. Yager, and S. Halloy for identifying the plants that grew in the field plots that received P. Funding: This work was supported by NSF grants for studying microbial community assembly in disturbed and periglacial environments (DEB-1258160 and PLR-1443578). Publication was funded by the University of Colorado Boulder Libraries Open Access Fund. Author contributions: S.K.S., C.C.C., D.R.N., J.E.K., and J.L.D. designed the experiments. J.L.D., S.K.S., J.E.K., S.C.C., C.C.C., and D.R.N. performed the experiments and wrote the paper. J.L.D. and S.K.S. performed bioinformatics and statistical analyses. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. DNA sequence data for this study have been deposited in the National Center for Biotechnology Information Sequence Read Archive, under BioProject ID PRJNA407491.

Stay Connected to Science Advances

Navigate This Article