Human disturbances alter the functioning and biodiversity of many ecosystems. These ecosystems may return to their pre-disturbance state after disturbance ceases, however humans have altered the environment in ways that may change the rate or direction of this recovery. For example, human activities have increased supplies of biologically-limiting nutrients, such as nitrogen (N) and phosphorus (P), which can reduce grassland diversity and increase productivity. We tracked the recovery of a grassland for two decades following an intensive agricultural disturbance under ambient and elevated nutrient conditions. Productivity returned to predisturbance levels quickly under ambient nutrient conditions, but nutrient addition slowed this recovery. In contrast, the effects of disturbance on diversity remained hidden for 15 years, at which point diversity began to increase in unfertilized plots. This work demonstrates that enrichment of terrestrial ecosystems by humans may alter the recovery of ecosystems and that disturbance effects may remain hidden for many years.

Experimental Design  

The work described here was conducted at the Cedar Creek Ecosystem Science Reserve (CDR) a U.S. Long Term Ecological Research (USLTER) site located (Latitude 45.4 N, Longitude 93.2 W) in Minnesota, USA (Figure S1). CDR lies on a sandy, outwash plain formed during the Wisconsin Glaciation about 11,000 years ago. These sandy soils are very low in N relative to other grasslands, and as a result have a relatively low productivity compared to other grasslands (Fay et al. 2015). The site has a mean annual precipitation of about 750 mm and a mean annual temperature of 6^o C (Borer et al. 2014a). During the course of this experiment, precipitation ranged from 538 to 1505 mm yr-1 (mean = 802 mm yr-1) and growing season (April – August) precipitation ranged from 281 to 970 mm yr-1 (mean = 483 mm yr-1). The experiment was replicated in three abandoned agricultural fields that were last tilled and farmed in 1968 (Field A), 1957 (Field B), and 1934 (Field C) (Figure S1). Detailed field descriptions are available in Tilman (1987).

The field experiment was established in 1982 and was composed of two treatments in a split-plot design: 1. Disturbance (Control or Disked, 35 x 55 m plots) and 2. Nutrient Addition (9 levels, 4 x 4 m plots). The nutrient-addition treatment had nine levels representing different combinations of Nitrogen (0 – 27.2 g N yr^-1 added as NH₄NO₃) and Other Nutrients (20 g m^-2 yr^-1 P₂0₅; 20 g m^-2 yr^-1 K₂0; 40 g m^-2 yr^-1 CaCO₃; 30.0 g m^-2 yr^-1 MgSO₄; 18 ug m^-2 yr^-1 CuSO₄; 37.7 ug m^-2 yr^-1 ZnSO₄; 15.3 ug m^-2 yr^-1 CoCO₂; 322 ug m^-2 yr^-1 MnCl₂; and 15.1 ug m^-2 yr^-1 NaMoO₄). Here we focus on two treatments: Control (No Nutrients added) and Fertilized (9.5 g N m^-2 yr^-1 added in combination with all other nutrients). Nutrients were applied twice per year in mid-May and mid-June. The complete list of treatments is presented in Table S1, and we present an analysis of all treatments in Figures S7-S10.

The Disturbance treatment was replicated in the three old-fields (A, B, C) in a completely randomized block design (2 treatments in each of 3 fields for a total of 6 35 x 55 m large plots). The nutrient treatments were replicated 6 times in a completely randomized design in each of the 35 x 55 m plots (54 4 x 4 m small plots). Thus, the complete experiment was composed of 3 fields, 6 35 x 55 m large plots, and 324 4 x 4 m small plots (Figures S1 and S2).

The establishment of the experiment proceeded as follows. Prior to the experiment, each of the fields was enclosed in a 1.8 m tall wire fence with 10 cm openings. In addition, woven wire fence with 6 mm openings was buried 84 cm in the ground and extended 60 cm above the ground. The fences reduced the density of mammalian herbivores, of which white-tailed deer (Odocoileus virginianus) and pocket gophers (Geomys bursarius) are of particular note at this site. Pocket gophers, in particular, can complicate interpretation of fertilization experiments, as they preferentially burrow in fertilized plots creating areas of bare soil and increasing spatial variability (Inouye et al. 1987; Inouye et al. 1997; Seabloom & Richards 2003). At this site, exclusion of deer and pocket gophers may have increased plant productivity, soil N, and colonization by woody plants (Inouye et al. 1987; Knops et al. 2000). It is unlikely that herbivore exclusion changed the relative effects of fertilization, based on experiments that manipulate both herbivore density and fertilization at our study site and other grasslands sites worldwide as effects tend to be additive (Inouye et al. 1987; Gruner et al. 2008; Borer et al. 2014b).

In April 1982, two 35 x 55 m areas were designated in each of the three fields (A, B, and C). In each of the fields, one of these two 35 x 55 m areas was selected to be disturbed with a 45 cm diameter disk harrow pulled by a tractor 20 times in one direction, 20 times perpendicularly, and 5 times diagonally to the first passes. Following the disking, the soil was hand raked to smooth the soil and remove any remaining vegetation, so that subsequent colonization was solely from seeds or small rhizome fragments. Within each of the 6 large plots, the 54 small plots were arrayed in 6 x 9 grid with 1 m buffers between each plot. Aluminum flashing was buried to depth of 30 cm around each plot to prevent horizontal movement of nutrients and spreading of plants through vegetative growth.

In our analyses, we focus on two nutrient treatments: 1. Control (no nutrients; Treatment I) and 2. Other Nutrients and 9.5 g of N (Treatment F)(Table S1). We chose the 9.5 g of N addition, as it is a rate that overcomes N limitation in our study systems and most grasslands systems without being toxic (Elser et al. 2007; Isbell et al. 2013a; Fay et al. 2015). We present analyses of the full N gradient in Figures S7 to S10 and Table S4. In brief, choosing a higher or lower N addition rate would change the strength of the N effect, but would not change the direction of the effects or change the strength of the interactions with the disking.

Sampling and Statistical Analyses

At peak biomass (mid-July to late August), all aboveground biomass was clipped in a 3 m by 10 cm strip (0.3 m²) in each plot. Note that there were 4 years when the disturbed plots were not sampled or only sampled in a single field. Inclusion or exclusion of years with missing data does not change the qualitative results. The biomass was sorted into dead, previous-year’s growth (litter) and live, current-year’s growth (live biomass). Live biomass was sorted to species, dried to constant mass at 40^o C, and weighed to the nearest 0.01 g. We estimated total aboveground biomass as the summed biomass of all non-woody species in each 0.3 m² sample, converted to g m^-2. We excluded woody biomass, because our goal was to estimate annual productivity and most of the woody biomass is from previous year’s growth. Woody plant biomass composed less than 1% of total biomass across the data set. Species richness is the number of species in each 0.3 m² sample. 

We quantified plant diversity as the Effective Number of Species based on the Probability of Interspecific Encounter (ENS_PIE), a measure of diversity that is more robust to the effects of sampling scale and less sensitive the presence of rare species than species richness (Jost 2006, 2007; Chase & Knight 2013). ENS_PIE is equivalent to the Inverse Simpson’s index of diversity which is calculated as where S is the total number of species (i.e., species richness) and p_i is the proportion of the community biomass represented by species i (Jost 2006, 2007; Chase & Knight 2013). Simpson’s evenness (E) satisfies the main requirements of an evenness index (Smith & Wilson 1996). In addition, it is directly related to ENS_PIE through the relationship E = ENS_PIE /S (Smith & Wilson 1996), thus we can factor diversity directly into its richness and evenness components through the relationship:

ENS_PIE = S * E

Across all data, ENS_PIE was positively correlated with richness (r=0.63) but uncorrelated with evenness (r=0.03). Richness and evenness were negatively correlated (r= -0.60).

All analyses were conducted in R (v. 3.5.3; R Foundation for Statistical Computing, Vienna, Austria). Mixed effects models were fit using the lmer function in the lme4 R library. We log₁₀ transformed all response variables prior to analysis to stabilize the variance of the residuals. In addition, the transformed data tested for proportional changes in the variables, which was biologically the most relevant given the among-field and among-year variability. In all models, among field and among year variability were treated as random effects. In the analyses of individual groups of species (e.g., annual or perennial plants), we analyzed log₁₀(biomass + 1), because there were some plots with zero mass of certain groups. All model specifications are included with each table of results. (Tables S5-S7).

Please refer to cdr-e001-e002-output-data-variable-descriptions for file descriptions.