Abstract Patterns of biodiversity are changing rapidly. “Legacy studies” use historical data to document changes between past and present communities, revealing long-term trends that can often be linked to particular drivers of ecological change. However, a single pair of historical samples cannot ascertain whether rates of change are consistent or whether the impact and identity of drivers have shifted. Using data from a second resurvey of 47 Wisconsin prairie remnants, we show that the pace of community change has increased with shifts in the strength of particular drivers. Annual rates of local colonization and extinction accelerated by 129 and 214%, respectively, between 1950 and 1987 and between 1987 and 2012. Two anthropogenic drivers—patch area and fire history—increased in importance between these periods. As the strength and number of anthropogenic forces increase, rates of biodiversity change are likely to accelerate in other ecosystems as well.

Keywords

Biodiversity

ecological communities

ecological diversity

ecology

environmental sciences

plant communities

novelty

INTRODUCTION As the world becomes increasingly dominated by humans (1, 2), ecosystems are being altered in unprecedented ways (3). Theory predicts that multiple global change drivers will lead to nonlinear rates of ecosystem change (4). However, these predictions are difficult to test empirically because of a lack of data, especially at large spatial and long temporal scales. One promising avenue for addressing these hypotheses is by means of “legacy studies” (5), in which historical data sets providing baseline records are revisited to characterize the magnitude of change over long periods of time (6–9). By substantially expanding temporal scales of ecological research, such legacy studies are capable of detecting long-term or lagged responses to drivers of change. However, many legacy studies are limited to a before-and-after comparison of two time points and, therefore, cannot assess theoretical predictions of nonlinearity. Understanding the pace of ecosystem change is motivated in part by the recognition that the number and strength of novel global change drivers are increasing (2, 4, 10). Although the distribution and dynamics of communities were previously known to be shaped by classic ecological drivers such as climate, soils, and species interactions, novel anthropogenic drivers such as habitat loss and disrupted disturbance regimes are increasingly interacting with or potentially outweighing these classic forces (4, 11). Legacy studies based on just two survey points not only prevent examination of the possibility of nonlinear changes but also prevent testing of whether indeed anthropogenic drivers are playing an increasingly important role in shaping ecological communities over time. We surmount these impediments by using a unique, multipoint data set that has the ability to detect accelerating rates of plant community change and the relative impact of environmental and anthropogenic drivers associated with this acceleration. These legacy data sets come from the University of Wisconsin Plant Ecology Laboratory (12) and include 47 remnant prairie sites surveyed in the 1950s [1947–1956 (13–15)], in the 1980s [1987–1988 (16, 17)], and again in 2012 (18). Temperate grasslands are productive and diverse ecosystems that are globally recognized for their contributions, such as carbon storage and agricultural production (19, 20). Worldwide, grasslands have disappeared more rapidly than any other biome (21), threatened primarily by conversion into agriculture (22) or shrubland (23). In North America, contemporary grassland remnants are largely confined to areas that are not suitable for agriculture, including rocky hill slopes, pioneer cemeteries, and railroad rights of way (13, 24). Disturbance by fire is an important driver of these prairie ecosystems, which are estimated to have burned historically every 1 to 3 years (25). Grazing by large ungulates has also been an important ecological and evolutionary force (26), though historical grazing in our study region is thought to be limited relative to other North American grasslands (27). In our study area, reduction in the frequency of fires began with European settlement in the 1830s; although small areas continued to be burned by private landowners or by sparks thrown by trains, overall fire frequency has declined dramatically (28, 29).

RESULTS AND DISCUSSION These prairie plant communities have considerably shifted in composition over the past six decades [permutational multivariate analysis of variance (PERMANOVA), F = 13.75, r2 = 0.17, P = 0.005; Fig. 1], distinguishing the contemporary communities from those that once occupied the same sites. Furthermore, the pace of change has accelerated. Plant community identity diverged more dramatically between 1987 and 2012 than between 1950 and 1987, despite the fact that the 1950–1987 interval is >30% longer. Annual rates of extinction increased by 214% between 1987 and 2012 relative to 1950–1987 (t = −10.11, df = 46, P < 0.001; Fig. 2), whereas annual rates of local colonization increased by 129% (t = −8.93, df = 46, P < 0.001). Despite these marked changes in species composition, the overall species richness has not changed significantly (F = 3.06, df = 2,138, P = 0.07). The average species richness of each site was 58, 59, and 51 in 1950, 1987, and 2012, respectively. Fig. 1 Species composition of prairie remnants has shifted over time. (A) Centroids for each survey are designated by a text box labeled with survey year. Blue ellipses are drawn at 1 SD from the centroid, and the dotted blue line describes the outer convex hull for each distribution. For full ordination, see fig. S1. The impact and identity of environmental drivers have also shifted. (B and C) The relationships between the three drivers and the changes in axis 1 (B) and axis 2 (C) scores are shown, with Spearman rank correlation (ρ) assessing the significance between changes in axis scores and environmental drivers. NMDS, nonmetric multidimensional scaling. Fig. 2 Colonizations and extinctions have increased. Annual rates of local colonization and extinction are both significantly higher between 1987 and 2012 than between 1950 and 1987 (paired t test, P < 0.001 in both cases). Error bars represent ±1 SE. The forces driving this acceleration in local extinction and colonization rates also appear to be changing. The influence of soil moisture (measured by the soil continuum index), a traditional environmental driver, has diminished relative to changes in fire history and patch area (both anthropogenic drivers). Between 1950 and 1987, the soil continuum index affected both colonization rate (F = 2.80, P = 0.038) and extinction rate (F = 2.69, P = 0.044; Fig. 3A), with fewer extinctions at drier sites. However, between 1987 and 2012, the soil continuum index was unrelated to local extinction rates, but local colonization rates were higher at wetter sites (F = 3.61, P = 0.013; Fig. 3B). In contrast, patch area was unrelated to annual extinction and colonization rates between 1950 and 1987 (Fig. 3C) but did affect annual extinction rates between 1987 and 2012 (r2 = −0.34, P = 0.02; Fig. 3D). As expected, the largest sites experienced the fewest extinctions. The changing role of patch area in this system matches the predictions of extinction debt (30, 31). Grasslands and other perennial plant communities often display time lags between the physical loss of habitat and ensuing species losses (11, 32, 33). Thus, although size was historically unrelated to species composition or dynamics, larger prairie patches now avoid more local extirpations and retain more of their historic species composition. Fig. 3 The impact and identity of the drivers of annual rates of colonization and extinction have changed. Between 1950 and 1987, colonization and extinction rates were significantly related to the soil continuum index, but not patch area (C) or fire history (F). Between 1987 and 2012, the soil continuum index (B), fire history (D), and patch area (F) were significant drivers of colonization and extinction. Trend lines in (C) to (F) represent statistically significant relationships. For statistical relationships, see table S1. Similar to patch area, fire history showed no relation to extinction and colonization rates between 1950 and 1987 (Fig. 3E). Fires were rare during this period, occurring at only 17% of sites. Since then, fires have become more common across all sites and were especially frequent in a subset of sites actively managed with prescribed fire. Sites with frequent fires experienced conspicuously fewer extinctions (r2 = −0.49, P < 0.001), more colonizations (r2 = 0.48, P < 0.001; Fig. 3F), and smaller shifts in floristic composition between 1987 and 2012. Other grassland studies suggested that fire maintains structure and diversity (16, 34, 35) by increasing habitat heterogeneity and by preventing expansion of woody species (36–38). Non-native species account for most of the colonizations that we documented at frequently burned sites between 1987 and 2012. This is an important and nonintuitive finding that merits further investigation, because frequent fire is generally thought to reduce the incidence and abundance of non-native species (39). One mechanism that may underpin this result is the reduction of accumulated litter following fire, which improves seed-soil contact and increases recruitment (40, 41). Both native and non-native species likely establish in the microsites opened by fire. However, our data document species occurrences rather than species cover. Thus, any additional recruitment of native species already present at a site will go undetected, whereas recruitment from a regional species pool increasingly dominated by non-native habitat generalists (42) would register as a colonization. Novel global change drivers may account for accelerated species turnover between 1987 and 2012 relative to 1950–1987. Patch size and fire frequency correspond to recent changes in these remnant prairie plant communities, whereas classic ecological drivers such as the soil continuum index are less related to community change. A similar shift toward anthropogenic global change drivers governing species composition has been documented in other systems (11, 43). The increased impacts of patch size and fire may be accompanied by changes brought about by other global change drivers. For example, long-term experiments demonstrate that grassland plant diversity generally declines in response to nitrogen deposition and elevated carbon dioxide (44–46). The Hierarchical Response Framework (HRF) predicts that these and other global change drivers resulting in long-term resource alterations may cause nonlinear ecosystem changes (4). Our findings demonstrate an accelerated pace of change consistent with the expectations of the HRF, using a different suite of drivers most relevant to our study system. Although further work is needed to disentangle the relative roles and potential interactions of numerous global change drivers, our results suggest that the pace of ecological change is increasing. The accelerated rates of change that we observe in these remnants reflect local extirpations and colonizations of particular groups of species. Most colonization events represent non-native species, which have increased in relative proportion across all sites by >500% between 1950 and 2012. At some sites, non-native species now account for >60% of total species richness. After the 1987 survey, Leach and Givnish (16) reported the steepest declines in species occurrence among nitrogen-fixing, short-statured, and small-seeded species, matching trends observed in other grassland studies (47–49); they attributed these losses to the lack of fire. We observed continuing declines in shorter plants and increases in taller species between 1987 and 2012 (fig. S2C and table S2). However, we did not observe continuing declines in smaller-seeded species or nitrogen-fixing species. Rather, the smallest-seeded species increased in frequency, whereas many larger-seeded species declined in frequency (fig. S2A and table S2). Although nitrogen-fixing species declined more rapidly than other species from 1950 to 1987, this trend disappeared between 1987 and 2012 (fig. S2B and table S2). Several factors might account for this, for example, increases in fires at some sites, declining rates of aerial nitrogen deposition, or the fact that losses of vulnerable nitrogen-fixing species had already played out by 1987. Rates of local colonization and extinction in these prairie remnants have accelerated by 129 and 214%, respectively, over the past six decades. Although the total species richness at each site remains similar, the floristic composition of species continues to diverge. The high rates of local extinction leave some sites with fewer than 18% of the species detected in the 1950 survey that are still present today. Meanwhile, we see a rapid increase in the rates at which non-native species are colonizing these sites, along with substantial shifts in the drivers associated with these changes. Local environmental factors that traditionally explained patterns of species composition (such as soil gradients) are declining in relative importance, whereas a suite of anthropogenic drivers are becoming more dominant both in these prairie remnants and in communities across the globe. Such altered relationships present difficulties for those entrusted to protect natural communities, as approaches and management strategies may be based on assumptions that no longer hold. However, we find some basis for optimism in the fact that some of the sites selected 60 years ago as the best remaining examples of Wisconsin prairies have largely retained their original plant community. This legacy study reveals that those prairies that most successfully avoid local extirpations and resist shifts away from historic species composition are large sites actively managed with regularly prescribed fire, providing empirical support for the assumptions that underpin conservation efforts in this system. Our results confirm the prediction that a nonlinear ecosystem response could be the new normal in an era of novel global change drivers. They also suggest that this response is best mitigated by restoring historic disturbances and by increasing available habitats.

MATERIALS AND METHODS Study system Our study sites are located in southern Wisconsin (geographic center, 43.23°N, 87.81°W), in the prairie-forest transition zone (fig. S3). The region includes both glaciated and unglaciated topographies, with elevations ranging from 177 to 488 m above sea level, and a continental climate. The mean annual temperature is 10°C, and the current mean annual precipitation ranges from 61 cm in the west to 84 cm in the east, which is an increase of 10 to 15% from 50 years ago (50). These changes have occurred regionwide, so we did not anticipate differential impacts of climate change among our study sites. Vegetation surveys Curtis and colleagues sampled >200 prairie remnants between 1947 and 1956 (referred to as the 1950 sample). Sites selected by Curtis (13, 14) and Curtis and Greene (15) for sampling had no evidence of recent human disturbance and, where possible, exceeded 1.6 ha (minimum size). For each site, presence/absence data were recorded by walking the entirety of the remnant. A subset of these remnants was sampled again in 1987 and 1988 (referred to as the 1987 sample) (16, 17). They relocated the sites using notes and maps from the original researchers and only sampled sites when there was little doubt that the correct location had been found, there was no sign of recent disturbance (that is, plowing or grazing by livestock), and the original prairie flora was reasonably intact. Fifty sites met the above criteria and were included in the 1987 sampling effort. Presence/absence data were recorded at each site by replicating the methodology of the original researchers. In the summer of 2012, we resampled the sites included in the 1987 survey (fig. S3). Three of the 50 sites sampled in 1987 were omitted because they could not be confidently relocated (one site) or had been severely disturbed (two sites; one was paved during road reconstruction and the other was subsumed into a homeowner’s lawn). Following the sampling methodology of previous researchers, we systematically walked the entire area of each site and recorded a list of all species present. We visited each site three times (in May, July, and August) to make a complete species list and to minimize the potential for intra-annual bias. Nomenclature was updated and standardized over the three sampling periods to account for changes in taxonomic classification. In instances where the original 1950 researchers did not identify plants to species level (for example, Carex and Rubus), we also considered these taxa at the genus level for analysis. We recognize that analyses of legacy data sets such as ours are sensitive to impacts from interannual variation. To quantitatively assess the magnitude of interannual variation in our system and to gauge the impacts on our findings, we resampled 10 of our sites in 2015 and compared these data to our original resampling effort in 2012. We did not find a substantial interannual variation in these remnants between 2012 and 2015 (table S3). Drivers of change To determine the relative impact of different types of drivers, we used a hypothesis-driven framework where we selected a priori drivers that were most relevant to our hypotheses and study system. We chose this approach, rather than a “kitchen sink” modeling framework where many potential drivers are evaluated, to prevent potentially spurious results. Using this framework, we examined the importance of changes in habitat patch area and fire history among time periods relative to the impacts of local soil moisture (a predominant classic driver of prairie plant community composition in our study area) (51) and more broadly (52, 53). We chose to evaluate changes in patch area because of the pervasive magnitude and known impact of habitat loss and fragmentation on the prairie biome in North America and within our study region (13, 22). In addition, we evaluated changes in fire regimes because of their documented importance for determining and altering prairie plant communities (16, 23). Although other global change drivers (for example, climate, nitrogen deposition, and increased CO 2 ) are known to alter the diversity and composition of grassland plant communities (44–46), they do not vary substantially among our study sites and thus were not included in our hypotheses or modeling framework. Patch areas were estimated in 1987 by pacing the perimeter of the patch (17) and in 2012 using field coordinates taken using a Global Positioning System unit with submeter accuracy (Trimble GeoXT). Because we lacked information on patch areas from the 1950 surveys, we quantified these from aerial photos taken between 1949 and 1956 and accessed from the Robinson Map Library at the University of Wisconsin. We georeferenced the photos in ArcGIS 10 (Esri, 2011) and calculated patch areas from polygons drawn around the extent of grasslands at each site. We corroborated the accuracy of the areas obtained through this methodology by comparing these to available data published for 23% of the sites (15) and found the estimates to be good approximations (r2 = 0.91). All patch areas were log-transformed to improve normality. We compiled fire histories for each site by interviewing local landowners, stewards, and managers following the methodologies of other studies (24, 54). Twelve of the 47 sites experienced prescribed fire in the past 60 years, ranging from 1 to 30 (mean, 9.25) fires during this time. Soil moisture was assigned using the soil continuum index devised by Curtis (54). Each site was given a position on a semiquantitative moisture scale based on the proportional abundance of select groups of indicator taxa (14, 15). The resulting five-unit ordinal index ranged from wet to dry, reflecting presumed soil moisture availability. A recent study using modern laboratory techniques examined several soil properties to evaluate the extent to which the Curtis index reflects reliable differences in soil moisture availability between sites (51). The study found that the index is indeed closely correlated with moisture supply, as measured by leaf δ13C. Data analysis All analyses were performed with R 3.1.2 (55). We used the vegan package (56) for multivariate procedures. Presence/absence data were used in all analyses. To determine whether community composition has changed among three surveys across 60 years, we used nonmetric multidimensional scaling (NMDS). This allows community composition changes to be visualized by creating a pair of orthogonal synthetic axes, on which sites can be plotted with proximity to other points indicating floristic similarity. To statistically evaluate significant differences in community composition across the three surveys, we used the PERMANOVA (57) function “adonis” in vegan, with Bray-Curtis distance as a measure of dissimilarity and 200 permutations. To assess whether rates of plant community change remained constant or were nonlinear, we compared the rates of local colonization and extinction over time within and among time periods. We calculated the annual rates of local colonization and extinction for each time period. We counted any species detected at a site that was not detected in the following resurvey as an extinction. Species undocumented in one survey but appearing in the subsequent resurvey were counted as colonizations. We used paired t tests to determine whether local colonization and extinction rates differed between 1950 and 1987 and between 1987 and 2012. We applied generalized linear models to assess the relationships between site characteristics (soil continuum index, patch area, and number of fires) and the annual colonization and extinction rates in each period (1950–1987 and 1987–2012). We log-transformed these response variables to meet model assumptions and matched the data on fire incidence and patch size to the time period (for example, colonization and extinction models from 1950 to 1987 used 1987 patch size data and counts of fire incidence for the 25 years preceding 1987, whereas the models for 1987–2012 used 2012 patch sizes and fire incidence in the 25 years preceding 2012). Values for the soil continuum index remained consistent in both models because we did not recalculate it, assuming that this site characteristic did not change between time periods (11). To understand the identity of the species driving the trends in community turnover, we looked at changes in occurrence based on seed size, plant height, and nitrogen-fixing ability. These traits are thought to be sensitive to fire suppression (16) and are reexamined here to facilitate comparison across time periods. We use trait data derived from the 1987 work and binned the continuous variables into the same classes used by Leach and Givnish (16) and Leach (17) (five classes for seed size and six classes for plant height). We tested the significance of proportional changes in functional identity using the McNemar change test with continuity correction (58)—a contingency test for paired data that is often used in re-survey studies [for example, in Leach and Givnish (16) and Keith et al. (59)].

SUPPLEMENTARY MATERIALS Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/2/2/e1500975/DC1 Fig. S1. Changes in species composition. Fig. S2. Effect of functional traits on species change. Fig. S3. Distribution of 47 prairie remnants in southern Wisconsin. Table S1. Results from four generalized linear models examining the role of site characteristics in annual colonization and extinction rates. Table S2. Changes in site occupation and net change among surveys as a function of plant traits. Table S3. Interannual variation in extinctions and colonizations since 1987 for 10 Wisconsin prairie remnants.

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 J. Curtis and all members of the Wisconsin Plant Ecology Laboratory for collecting and archiving the historical data. Funding: Funding was provided by Prairie Biotic Research Inc., the University of Wisconsin–Madison Zoology Department, and the Kettle Moraine Garden Club Scholarship Fund. Author contributions: The 1987 data were collected by M.K.L. and J.A.H. The 2012 data were collected by A.O.A. The Damschen Laboratory provided assistance with data collection and writing of the paper. Authorship was determined by the sequence-determines-credit approach with alphabetical listing, where equal contributions were made. A.O.A. and E.I.D. led, and all other authors contributed to idea generation, analysis, and writing of the paper. 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 and can be found in Dryad:10.5061/dryad.jc00c.