Significance Human activities have elevated nitrogen (N) deposition and there is evidence that deposition impacts species diversity, but spatially extensive and context-specific estimates of N loads at which species losses begin remain elusive. Across a wide range of climates, soil conditions, and vegetation types in the United States, we found that 24% of >15,000 sites were susceptible to N deposition-induced species loss. Grasslands, shrublands, and woodlands were susceptible to species losses at lower loads of N deposition than forests, and susceptibility to species losses increased in acidic soils. These findings are pertinent to the protection of biodiversity and human welfare and should be considered when establishing air quality standards.

Abstract Atmospheric nitrogen (N) deposition has been shown to decrease plant species richness along regional deposition gradients in Europe and in experimental manipulations. However, the general response of species richness to N deposition across different vegetation types, soil conditions, and climates remains largely unknown even though responses may be contingent on these environmental factors. We assessed the effect of N deposition on herbaceous richness for 15,136 forest, woodland, shrubland, and grassland sites across the continental United States, to address how edaphic and climatic conditions altered vulnerability to this stressor. In our dataset, with N deposition ranging from 1 to 19 kg N⋅ha−1⋅y−1, we found a unimodal relationship; richness increased at low deposition levels and decreased above 8.7 and 13.4 kg N⋅ha−1⋅y−1 in open and closed-canopy vegetation, respectively. N deposition exceeded critical loads for loss of plant species richness in 24% of 15,136 sites examined nationwide. There were negative relationships between species richness and N deposition in 36% of 44 community gradients. Vulnerability to N deposition was consistently higher in more acidic soils whereas the moderating roles of temperature and precipitation varied across scales. We demonstrate here that negative relationships between N deposition and species richness are common, albeit not universal, and that fine-scale processes can moderate vegetation responses to N deposition. Our results highlight the importance of contingent factors when estimating ecosystem vulnerability to N deposition and suggest that N deposition is affecting species richness in forested and nonforested systems across much of the continental United States.

Global emissions of reactive nitrogen (N) to the atmosphere and subsequent deposition into terrestrial ecosystems have tripled in the last century (1). This N deposition has been identified as a threat to plant diversity (2⇓–4), and plant diversity is linked to ecosystem stability (5), productivity (6), and other ecosystem services (7). Elevated nitrogen inputs have been shown to cause decreases in species richness over time in small plot experiments (8⇓–10) and in regional gradient studies in Europe (11, 12). Although these studies and others have led to some generalizations about the impacts of N deposition on plant diversity, most of these studies have focused on grassland ecosystems and/or, in the United States, have been fine-scale field experiments where N is added experimentally as fertilizer. Thus, translation of these findings to nongrassland systems or to large regions of the United States may not be appropriate. Unlike grasslands, where elevated N has often led to light limitations and subsequent competitive exclusion (13), plant growth in the herbaceous layers of forest understories is typically primarily light-limited (14) regardless of the extent of N inputs. Moreover, soil chemistry can be heterogeneous, influencing the potential of soil acidification by nitrogen deposition (15). In most arid ecosystems, moisture may be more important than nutrients in controlling plant growth during the growing season (16, 17). Finally, the level of N input at which diversity is first impacted (18) is often unknown for many regions because most studies use a fairly coarse experimental approach to estimate thresholds of response or have been conducted where there have already been high inputs of N for decades (e.g., Northern Europe). To address these critical gaps in our knowledge of continental-scale relationships between N deposition and plant diversity, we used data from herbaceous ground-layer communities within 15,136 forest, woodland, shrubland, and grassland sites spanning N deposition gradients across the continental United States. More specifically, we assessed how covarying climate and edaphic factors affected ecosystem vulnerability to N deposition.

Nitrogen inputs can increase diversity, decrease diversity, or leave diversity unchanged, contingent on a host of associated ecosystem factors. Biodiversity can be reduced through several general mechanisms (4), including but not limited to (i) release from N limitation that leads to increased aboveground production, reduced light availability, and ultimately competitive exclusion (13, 19) and (ii) soil acidification and associated cation depletion and imbalances that lead to recruitment inhibitions (20, 21). The importance of N limitation likely declines in arid areas that are more moisture-limited or in warm, wet areas favoring high net N mineralization, either one of which may reduce the importance of external N inputs. In such cases, N may be less limiting to plant growth, and therefore communities are less responsive to additional N deposition (2). Conversely, enrichment may increase biodiversity in extremely N-poor environments where release from N limitation does not result in competitive exclusion (22, 23) or where soils have a high pH resistant to soil acidification (11, 24).

Because N enrichment can affect plant diversity through multiple pathways and environmental contingencies, we investigated whether N deposition is a widespread threat to plant species diversity or whether some vegetation types or environments are more vulnerable than others. We compiled herbaceous plant species composition data from existing datasets (Table S1) that included 15,136 sites and 3,852 herbaceous species from across the continental United States. At each site, we calculated species richness, the total number of unique species per plot, a commonly used metric of diversity (25). We then extracted geospatial estimates (Table S2) of N deposition, annual precipitation, mean annual temperature, and soil pH for each site. As in several previous studies in Europe (11, 12, 26), we used a correlative approach that cannot show direct causality but can nevertheless provide insight into the mechanisms involved in, and communities most susceptible to, loss of diversity as a result of N deposition. First, we analyzed relationships between plant species richness and N deposition involving interactions with precipitation, temperature, and soil pH within two broadly defined vegetation types (closed canopy forest vs. open canopy grasslands, shrublands, and woodlands). We then examined the same set of predictors within gradients defined by unique combinations of specific vegetation communities and source datasets that spanned an adequate range of N deposition (Methods and Table S3).

Table S1. Summary of vegetation data sources for national surface analysis

Table S2. Summary of predictor variable sources

Table S3. Vegetation Alliances with sufficient data for gradient analyses

Conclusion Our continental-scale analysis found that the threat of N deposition to herbaceous plant species richness is ecosystem-specific, with some ecosystems more vulnerable than others, and some conditions conferring greater vulnerability. Ecosystems with open vegetation (grasslands, shrublands, and woodlands) had lower critical loads of N deposition (7.4–10.3 kg N⋅ha−1⋅y−1) than ecosystems with closed-canopy forest vegetation (7.9–19.6 kg N⋅ha−1⋅y−1). Within these broad vegetation groups, declines in species richness along gradients of increasing N deposition were more likely to occur in ecosystems with acidic soils. Climate also interacted with N deposition to help explain species richness, but its influence was less consistent across scales. Increasing the number of N-addition experiments with treatment levels spanning 2–20 kg⋅ha−1⋅y−1 and implementing them across the full range of soil pH, climate, and vegetation types that exist on the landscape would be a very welcome complement to the correlative work that we have reported here. In the meantime, our work suggests that the mechanism of competitive exclusion via shading is likely of reduced strength in the comparative shade of forest understories whereas the acidification and competitive exclusion mechanisms are probably more likely to occur synergistically in the high-light environment characteristic of grasslands. We successfully identified ecosystems vulnerable to N deposition and refined herb-based N deposition critical loads (18) by incorporating a broad range of vegetation types, N deposition loads, soil substrates, and climate conditions in our analysis. This identification of vulnerable ecosystems and influential environmental factors is critical for managers to set monitoring and conservation priorities.

Methods Data Acquisition and Management. We compiled vegetation data from multiple sources (Table S1) because a single standardized national dataset of herbaceous plant species presence and abundance with sufficient spatial coverage and plot density is not available for the United States. We retained only terrestrial sites sampled after 1989 that had a complete inventory of species from graminoid and forb functional groups, quantitative abundance for each plant species, a sampling area of 100–700 m2, and known geographic coordinates. At each site, we calculated total herbaceous (defined here as forbs and graminoids) plant species richness, a conservative measure because total richness could remain unchanged even as invasive species richness increases and native species richness declines. We estimated N deposition by adding Community Multiscale Air Quality (CMAQ) model dry deposition estimates to interpolated National Atmospheric Deposition Program (NADP) wet deposition and extracting a value based on coordinates for each site. The CMAQ version 5.0.2 dry deposition estimate was a 10-y average (2002–2011) with 12-km resolution, using models run in 2014 by Robin Dennis at the Environmental Protection Agency (EPA). CMAQ dry deposition estimates, or other comparable estimates with fine resolution, are not yet available at a national scale before 2002. The NADP wet deposition was a 27-y average (1985–2011), which we resampled from the raw 2.33833-km resolution to the 4-km resolution of the Parameter-Elevation Relationships on Independent Slopes Model (PRISM) precipitation data that had been used in the interpolation. We extracted climate covariates [specifically, average annual precipitation and temperature from 30-y PRISM climate normals (1981–2010)] and obtained soil pH, where available, from the same datasets that supplied vegetation data. If soil data from soil samples colocated with vegetation data were not available, then pH from 1:1 water extracts from the national US Department of Agriculture (USDA) Soil Survey Geographic (SSURGO) database was used. We retained the 15,136 sites with nonmissing species richness and predictor values that met the criteria for analyses at either the national scale (data sources combined but plots filtered based on area) or gradient scale (data sources considered separately). Data Analysis. For our initial national-scale analysis, we began with all 15,136 sites, and then, based on expected differences in mechanisms, we divided those sites into two broad vegetation types: namely, closed canopy (deciduous forest, evergreen forest, and mixed forest) and open canopy (grassland, shrubland, and woodland) vegetation types. Within each of these two groups, we determined the relative importance of our four primary predictor variables (N deposition, soil pH, precipitation, and temperature) by looking at the R1 coefficients of determination (based on absolute deviations in quantile regression rather than squared deviations) of b-spline models with and without these four main effects. Next, we examined nonlinear regressions of the 0.50 (median), 0.10, and 0.90 quantiles of total herbaceous plant species richness response to N deposition (quadratic), soil pH, mean annual temperature, annual precipitation, and the two-way interactions involving N deposition (i.e., N × precipitation, N × temperature, and N × soil pH) using the quantreg package of R (version 3.0.2) software. Out of all possible models, we selected the model with the lowest corrected Akaike information criterion (AICc) for each of the two broad vegetation types (Table 1 and Fig. 1). We used the median quantile regression model with the best AICc to calculate separate critical loads of N deposition for open and closed canopy vegetation. Qualitatively, critical loads of N deposition are defined here as the N deposition threshold at which species richness begins to decline, corresponding graphically to the N deposition level at which a hump-shaped relationship between N deposition and species richness reaches its peak value of species richness. Quantitatively, we calculated critical loads of N deposition by taking the first derivative of the best model with respect to nitrogen and setting that expression to zero, for models with a negative quadratic N deposition term. For critical loads specific to each site, we used the coefficients from the critical load expression and site-specific covariate values. We subtracted critical loads from N deposition to determine exceedances of N deposition critical loads. Three sets of exceedances were calculated, using (i) the median point estimates of critical loads, as well as (ii) the upper and (iii) the lower limits of the 95% CI of the critical loads. Only the exceedances based on the median point estimates of critical loads are presented graphically and in the Abstract. Further community-scale analyses were focused on individual alliances as defined by the National Vegetation Classification (NVC) (33). We analyzed alliances with deposition gradients with maximum N deposition that was either 2.5 times or 4 kg⋅ha−1⋅y−1 greater than minimum N deposition, and that had at least 20 sites from at least one common data source. These gradient criteria reduced the number of sites to 6,807. For each N deposition gradient, we performed multiple regressions of species richness against N deposition, with the same predictor variables and the same model selection procedure as in the national analysis (except that N deposition was only first order).

SI Methods Data Acquisition and Management. Vegetation data sources included the US Forest Service (USFS) Forest Inventory and Analysis Phase 3 Vegetation Indicators herbaceous plant dataset, the Ecological Society of America’s VegBank archive of vegetation plots (www.vegbank.org), state Natural Heritage Programs, individual researchers, and other organizations (Table S1). Vegetation datasets that were themselves collections of datasets (MN and VA) were subdivided based on the agency or project that collected the data. Species-area curves were essentially flat, presumably because sites had already mostly saturated by 100 m2. We corrected misspelled plant species names, revised taxonomy to match use at USDA PLANTS (plants.usda.gov/java/), lumped subspecies and varieties to the species level, and added plant species attributes obtained from USDA PLANTS (plants.usda.gov/adv_search.html). Each vegetation sampling site was classified into a vegetation type: specifically, the alliance level of the 1997 version of the National Vegetation Classification (NVC) (33), which generally specifies the one to three most dominant plant species. Where investigators had used a classification scheme other than NVC alliances, we reassigned sites to the closest possible alliance. Unclassified unforested sites were assigned an alliance value based on dominant species within the field vegetation data. Our modern (1985–2011) N deposition estimate correlated well with a short 5-y temporal subset of N deposition from 2006 to 2010 [r2 = 0.98, Ndep 2006–2010 = (Ndep 1985–2011 *0.84) + 0.53], with the 2006–2010 subset averaging 0.80 kg⋅ha−1⋅y−1 lower. Critical loads of N deposition based on just 5 y of N deposition likely underestimate the importance of N retention over multiple decades whereas critical loads based on N deposition accumulated over the course of 160 y would likely underestimate the importance of N losses. In contrast, the intermediate timescale of our N deposition estimate (10 y of dry deposition and 27 y of wet deposition) is consistent with a moderate and realistic amount of N retention. Nevertheless, to illustrate the possible role of temporal trends in N deposition, we conducted a parallel analysis using just the 5-y temporal subset of N deposition from 2006 to 2010, finding that exceedances of N deposition critical loads calculated using the 2006–2010 N deposition subset were on average just 0.83 kg⋅ha−1⋅y−1 higher than exceedances calculated using the 1985–2011 N deposition estimates [r2 = 0.97, Exceedance 2006–2010 = (Exceedance 1985–2011 *0.82) + 0.11]. Additional alternative N deposition estimates were considered. Our modern (1985–2011) N deposition estimate correlated well (r2 = 0.78) with average N deposition that also includes historic estimates from 1850 to 1984 (34) [Ndep historic = (0.63*Ndep 1985–2011 ) − 0.02], despite the coarse scale (∼2° latitude × 5° longitude grid) and uncertainty of the historical data. Total CMAQ N deposition was correlated with our total N deposition estimates [r2 = 0.98, CMAQ 2002–2011 = (Ndep 1985–2011 *0.99) − 0.24]. An initial version of total deposition estimates from the Total Deposition Science Committee (TDEP) (35) was correlated with our total N deposition estimates [r2 = 0.89, TDEP 2000–2012 = (Ndep 1985–2011 *0.91) + 0.30], and, when TDEP estimates are finalized, they will be useful for some future analyses. Monitoring data and spatial modeling techniques are not yet sufficiently reliable to add organic nitrogen to our inorganic nitrogen deposition estimate. We also examined housing density as a potentially confounding indicator of urbanization but found that its correlation coefficient with N deposition was only 0.19 and did not use it in further analyses. The housing data were those data used in Radeloff et al. (36). Data Analysis. For the national surface analyses, we characterized uncertainty in the site-specific critical load values by performing a Monte Carlo resampling of the coefficients from a multivariate normal distribution, with means equal to the coefficient estimates and variance/covariance from the parameters in the median quantile regression model. We computed a 95% confidence interval for each site estimate of critical loads based on the 2.5th and 97.5th percentiles of 10,000 simulations for each site. For the gradient analyses, we extracted the model coefficients (and their 95% CI), the r2, and the minimum and maximum values of each predictor variable for each gradient. Within each gradient, we calculated species richness change as the simple slope of N deposition from multiple regression coefficients: β N + (β N*M × M i ), where β N is the parameter for N deposition, β N*M is the parameter for the interaction of N deposition and the moderating variable M, and M i are values of the moderating variable M across the gradient. Using the N deposition slopes calculated above and the observed range of moderating variables in each gradient, we classified each gradient as containing one of the following relationships between N deposition and plant species richness: (i) an unconditionally negative relationship, (ii) an unconditionally positive relationship, (iii) one of several specified relationships that are conditional on other moderating environmental predictors, or (iv) no detectable relationship between N deposition and plant species richness. For gradients that had a detectable relationship between N deposition and plant species richness, we also plotted the N deposition simple slopes as a function of the moderating variables of each gradient.

Acknowledgments Vegetation data were shared by the Forest Inventory and Analysis Database (FIADB) Vegetation Indicators Program, the Ecological Society of America VegBank, the Minnesota Biological Survey, the New York, Virginia, and West Virginia Natural Heritage Programs, Robert Peet and the Carolina Vegetation Survey, the US National Park Service Southern Colorado Plateau Network, the University of Wisconsin Plant Ecology Laboratory, Kevin Knutson of the US Geological Survey (USGS), and the coauthors. This paper arose from the “Diversity and Nitrogen Deposition” working group supported by the John Wesley Powell Center for Analysis and Synthesis, funded by the USGS. The US Environmental Protection Agency (Contract EP-12-H-000491) and the Cooperative Ecosystem Studies Units Network (National Park Service Grant P13AC00407 and USGS Grant G14AC00028) provided additional funding. The USGS supports the conclusions of research conducted by their employees, and peer reviews and approves all of their products consistent with USGS Fundamental Science Practices. The views expressed in this manuscript do not necessarily reflect the views or policies of the US Environmental Protection Agency or the USDA Forest Service. Mention of trade names or commercial products does not constitute endorsement or recommendation for use.

Footnotes Author contributions: S.M.S., E.B.A., W.D.B., C.M.C., J.B., and M.L.B. designed research; S.M.S., E.B.A., W.D.B., and C.M.C. performed research; S.M.S. and C.M.C. analyzed data; and S.M.S., E.B.A., W.D.B., C.M.C., J.B., M.L.B., B.S.C., S.L.C., L.H.G., F.S.G., S.E.J., L.H.P., B.K.S., C.J.S., K.N.S., H.L.T., and D.M.W. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

Data deposition: The data reported in this article have been deposited in the Dryad Digital Repository, datadryad.org (doi: 10.5061/dryad.7kn53).

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1515241113/-/DCSupplemental.