Examining the assumptions of heterogeneity-based management for promoting plant diversity in a disturbance-prone ecosystem

Study site

We conducted our study at the Tallgrass Prairie Preserve (TGPP) which is a 15,700 ha nature preserve that is managed using a patch-burn approach (Hamilton, 1996, 2007; Allen et al., 2009). The TGPP is located between 36.73° and 36.90° N latitude, and 96.32° and 96.49° W longitude, in Osage County, Oklahoma and owned by The Nature Conservancy (TNC). Over the course of the 12 year study period (1998–2009), total annual rainfall varied from 494 to 1,252 mm. The preserve is situated at the southern extent of the Flint Hills region. The elevation of the preserve ranges from 253 to 366 m, and the underlying bedrock of the region is characterized by soils derived from Permian limestone, sandstone, and shale (Oviatt, 1998). Due to the proximity of bedrock to the surface and the relatively steep terrain, the Flint Hills region has experienced long-term erosion leaving surface layers of soil that are thin and young. Because of this rockiness the Flint Hills region, including the TGPP, has remained unplowed and has been instead used primarily as rangeland for cattle. Prior to the acquisition of the preserve by TNC in 1989, the majority of the site was managed for cow-calf and yearling cattle production with a 4- to 5-year rotation of prescribed burning and patchy aerial application of broadleaf herbicides (1950–1989) (Hamilton, 2007).

Approximately 90% of the TGPP consists of grasslands. The majority of the grasslands are composed of tallgrass prairie habitats dominated by Andropogon gerardii, Sorghastrum nutans, Sporobolus compositus, Panicum virgatum, and Schizachyrium scoparium. Shortgrass prairie habitat occurs to a lesser extent on more xeric sites and is dominated by Bouteloua spp. Despite the application of herbicide earlier in the 20th century, the flora of the preserve appears relatively intact with a total of 763 species of vascular plants (to date) of which 12.1% are exotic (Palmer, 2007).

Disturbance regime

The management plan at the TGPP encompasses a wide range of spatial and temporal variation in the application of prescribed fire and cattle or bison grazing (Hamilton, 1996, 2007). In 1993, 300 bison were introduced year-round onto a 1,960 ha portion of the preserve. As the bison herd increased in size, the area allotted to the herd was increased eight times to an area of 8,517 ha by 2007 (Fig. 1A, 54% of preserve area). Initial bison stocking rates were increased in 1999 to 2.1 animal-unit months ha⁻¹ (see Hamilton, 2007 for additional details). Within the bison unit, animals were allowed to range freely and their movement was not obstructed by internal fences. Watersheds within the bison unit were considered randomly for burning only if they met the minimum fuel criteria of 900 kg ha⁻¹ of fine fuels. Within a given year, the season of burn of the bison unit was split as follows: 40% dormant spring (March–April), 20% late growing season (August–September), and 40% dormant winter (October–December). The remainder of the preserve was seasonally grazed by cattle and typically burned more frequently in the dormant spring season, but some of the cattle pastures were utilized for smaller scale (2,350 ha) patch-burn experiments in which only one-third of a given management unit was burned annually (Hamilton, 2007). Stocking within the cattle pastures included both intensive-early stocking and season-long stocking, which contrasted with the year-round stocking in the bison unit.

Data collection

During a 3-year period (1997–2000) we collected vegetation and environmental data at every intersection of a one km UTM grid for a total of 151 samples (Fig. 1B). For our analysis we only used 128 of the samples which had no standing water and less than 20% of cover by rocks or woody plants to ensure that they were representative of grasslands. We randomly selected 20 of these plots for annual resampling, which continued from 1998 to 2009 (Fig. 1C). All vegetation samples were collected in June (when we could readily identify both early and late-season plants). We used a plot design modified from Peet, Wentworth & White (1998). At each site we recorded the presence/absence for plants identified to the species level that were rooted within each of the four corners of the 10 × 10 m quadrats at four nested spatial scales: 0.01, 0.1, 1, and 10 m² (see McGlinn, Earls & Palmer, 2010 for a sampling diagram). At the 100 m² spatial scale we visually estimated species percent cover using cover classes (1: trace, 2: <1%, 3: 1–2%, 4: 2–5%, 5: 5–10%, 6: 10–25%, 7: 25–50%, 8: 50–75%, 9: 75–100%) rather than continuous estimating cover to decrease the amount of error in our cover estimates (Peet, Wentworth & White, 1998). For the analysis, we used the mid-point of each cover class. For the purposes of this study we only report results at the 100 m² scale, but analyses at the finer spatial grains (using the corner subplots) yielded qualitatively similar results but higher levels of unexplained variance.

At every sampling event we combined four 15 cm soil cores collected at each corner of the plot. We sent the soil samples to Brookside Labs (New Knoxville, OH, USA) to be analyzed for soil cations: P, Ca, Mg, K, Na, B, Fe, Mn, Cu, Zn, and Al. The resampled vegetation and environmental data are in a public online archive (McGlinn, Earls & Palmer, 2010).

Data analysis

Our analysis was composed of two sets of parallel analyses: those on species richness and those on species composition. Additionally, we carried out separate analyses on the 128 vegetation samples which were sampled once over a 3-year period and the 20 samples that were annually resampled over a 12-year period (n = 20 × 12 = 240). We will refer to these separate analyses as the Grid analysis and the Repeat analysis, respectively. The Grid analysis was used to test assumption I and the Repeat analysis was used to test assumption II of the patch-burn approach.

Assumption I asserts that spatial variation in the community should reflect spatial variation in management. To carry out a strong test of this hypothesis we compared the effect sizes and variance explained by management variables relative to inherent environmental heterogeneity captured by soil cations on plant community richness and composition. We quantified variation in soil using the first three axes of a principal components analysis on the 11 soil cations (see Table S1 for cation PC loadings). We used three axes because these explained the majority of the variation in the dataset (68%) and this provided a balanced comparison of three soil variables against three management variables. We chose to examine soil variables because previous research at the site indicated that they are good proxies for spatially relevant environmental variation (Palmer et al., 2003).

Assumption II asserts that local patches change through time in response to disturbance. To test assumption II we examined the degree to which spatio-temporal changes in richness and composition correlated with changes in management after controlling for inherent site and year specific differences. While the Grid analysis examined the effects of key soil variables on species composition, it is possible that unmeasured variables that distinguish one site from another are responsible for intersite differences. Similarly, interannual variation may be inadequately summarized by readily obtainable climatic variables. Thus, the simplest way to completely account for interplot and interyear differences, independent of management, is to use site and year nomial variables, respectively.

In all the analyses, we quantified management using three variables: years of bison management, years since burn, and number of burns in the past 5 years. Year of bison management is meant to capture variation attributed to the differences in bison and cattle management at the preserve. Note that sites not grazed by bison were grazed by cattle and there were no ungrazed sites. We chose not to include season of burn as an explanatory variable because 83% (67 out of 80) of the prescribed fire events recorded on our study sites took place during the dormant season.

In both the Grid and Repeat analyses, we quantified the variance explained by the competing classes of variables using variation partitioning which estimates the unique and shared fractions of explained variance in species richness and composition (Peres-Neto et al., 2006; Legendre & Legendre, 2012). The univariate species richness analyses were carried out using ordinary least squares regression (OLS). We also carried all of our analyses out using generalized least squares (GLS) regression that accounted for spatial and temporal autocorrelation in the error of our models (Pinheiro & Bates, 2000), but we found that the GLS model had qualitatively similar results to the OLS model. Therefore, we only report results from the OLS analysis below.

We analyzed species composition using redundancy analysis (RDA) and canonical correspondence analysis (CCA) direct ordination approaches. We only report the results of the RDA analyses because they were qualitatively similar to the CCA results. For the compositional analyses, we down-weighted the effect of common species using a square root transformation of the cover estimates. Following the recommendations of Peres-Neto et al. (2006), we report the adjusted coefficient of determination R_adj² for each fraction of variance using Ezekiel’s adjustment Ezekiel (1930) for the OLS and RDA analyses. We tested if the individual fractions of variation were significantly different than zero using permutation tests with 999 permutations of the rows of the response variable of interest.

All analyses were conducted in R (R Development Core Team, 2018) and the multivariate analyses were carried out using the package vegan (Oksanen et al., 2018). The code and data are publically available at an online repository (https://zenodo.org/record/2641111#.XLUvR-hKhnI). The data shared at this repository includes the data shared by McGlinn, Earls & Palmer (2010), two-additional years of data, and the vegetation, management, and environmental data for the remainder of the plots on the UTM grid.

Grid analysis—test of spatial management effects

The OLS model explained 17% of the total variance in richness (F_6,121 = 4.13, p < 0.001), and the RDA model explained 16% of the total variance in species composition (F_6,121 = 3.76, p < 0.001). Spatial variation due to differences in soil properties explained the largest proportion of variance in species richness and composition (Fig. 2, R_OLSadj² = 0.15 and Fig. 3, R_RDAadj² = 0.08, respectively). The first soil PC axis reflected a gradient from high iron sites to high calcium sites (Table S1), and it was the most strongly correlated variable with richness (standardized partial regression coefficient, β = −0.38, p < 0.001, Table S2) and composition (Table S3). The independent effect of management was negligible on species richness (Fig. 2, R_OLSadj² ≊ 0, p = 0.59) but not on species composition (Fig. 3, R_RDAadj² = 0.02, p = 0.001).

Repeat analysis—test of space-time management effects

For the repeat analysis, the OLS model explained 76% of variance in richness (F_33,206 = 19.23, p < 0.001), and the RDA model explained 61% of the total variance in species composition (F_33,206 = 9.78, p < 0.001). We observed statistically significant management effects on both richness (Fig. 4) and composition (Fig. 5); however, relative to site and year effects management explained little of the variance (R_OLSadj² = 0.04 and R_RDAadj² = 0.01 for richness and composition, respectively, Fig. 6). The majority of variance in richness and composition was attributed to site specific effects (R_OLSadj² = 0.49 and R_RDAadj² = 0.45, respectively, Fig. 6). Year effects were intermediate to those of site and management effects. The majority of the variance attributed to management was shared with site effects (Fig. 6). Years of bison management was the most important management variable (Figs. 4 and 5; Tables S4 and S5) which was positively correlated with richness (β = 0.46, p < 0.001, Fig. 4). Years of bison management increased cover of many weedy annual plants such as Amphiachyris dracunculoides (amphidrac), Kummerowia stipulacea (kummstip), Chamaecrista fasciculata (chamfasc), Medicago lupulina (medilupu), and Melilotus officinalis (melioffi). Bison management was correlated with decreased cover of several long-lived perennials species such as Helianthus mollis (helimoll), Asclepias viridis (asclviri), Spiranthes spp. (unkspir), and Viola sororia (violsoro) (Fig. 5). The number of burns in the past 5 years and time since the last burn were both significantly negatively correlated with species richness (Fig. 4). Frequent burning increased Schizachyrium scoparium (schizscop) but decreased Andropogon gerardii (andrgera) and Sporobolus compositus (sporcomp) (Fig. 5).