Analysis of abrupt transitions in ecological systems

. The occurrence and causes of abrupt transitions, thresholds, or regime shifts between ecosystem states are of great concern and the likelihood of such transitions is increasing for many ecological systems. General understanding of abrupt transitions has been advanced by theory, but hindered by the lack of a common, accessible, and data-driven approach to characterizing them. We apply such an approach to 30–60 years of data on environmental drivers, biological responses, and associated evidence from pelagic ocean, coastal benthic, polar marine, and semi-arid grassland ecosystems. Our analyses revealed one case in which the response (krill abundance) linearly tracked abrupt changes in the driver (Pacific Decadal Oscillation), but abrupt transitions detected in the three other cases (sea cucumber abundance, penguin abundance, and black grama grass production) exhibited hysteretic relationships with drivers (wave intensity, sea-ice duration, and amounts of monsoonal rainfall, respectively) through a variety of response mechanisms. The use of a common approach across these case studies illustrates that: the utility of leading indicators is often limited and can depend on the abruptness of a transition relative to the lifespan of responsive organisms and observation intervals; information on spatiotemporal context is useful for comparing transitions; and ancillary information from associated experiments and observations aids interpretation of response-driver relationships. The understanding of abrupt transitions offered by this approach provides information that can be used to manage state changes and underscores the utility of long-term observations in multiple sentinel sites across a variety of ecosystems.


INTRODUCTION
Many ecological systems can exist in two or more states that differ in abundance or composition of species, rates of ecological processes, and ecosystem services provided by them (Beisner et al. 2003, Suding et al. 2004. Smooth, gradual transitions between ecosystem states are unremarkable, occurring during succession or as ecosystems track gradually changing environmental conditions. In contrast, abrupt transitions between ecosystem states are typically unexpected and can have wide-ranging, negative impacts. Abrupt transitions happen either when the gradually changing environment passes a critical point or when discrete perturbations cause sudden changes in underlying environmental drivers. Abrupt and irreversible transitions are forecast to increase as climatic changes and depletion of natural resources both accelerate (Millennium Ecosystem Assessment 2005, Fagre et al. 2009). Such forecasting, however, is difficult because there are many different causes of state changes (Hastings and Wysham 2010) and because existing approaches demand far more data than are normally available (Carpenter et al. 2011).
Managing state changes is as difficult as forecasting them. When environmental changes are not severe, or when organisms with short lifespans and generation times rapidly track environmental drivers, some state changes can be reversed in relatively short periods of time ( 50 years) if drivers are returned to pre-change conditions or perturbations are eliminated (Jones and Schmitz 2009). In other cases, environmental change can result in state changes that persist long after environmental drivers have returned to earlier conditions. The persistence of these socalled ''ecological thresholds'', ''regime shifts'', ''phase shifts'', or ''catastrophes'' (Hughes 1994, Scheffer et al. 2001, Groffman et al. 2006) is caused by time-lags in the responses of biological systems to environmental change (hysteresis), differences in the relationships between state variables and environmental drivers before and after the state change, or the appearance of novel feedbacks among state variables and drivers that reinforce the new state (Scheffer et al. 2001, Lindig-Cisneros et al. 2003, Briske et al. 2006, Suding and Hobbs 2009).
The development of management strategies to mitigate abrupt transitions requires strong linkages among theory, data, and case studies, but there is little guidance available for using historical or ongoing studies to detect or respond to abrupt transitions. There is confusion and disagreement about what changes constitute transitions Davis 2003, Schroder et al. 2005) and a limited understanding of ecological mechanisms causing them Archer 1999, Collie et al. 2004). Empiricists disagree about how to best gather and interpret relevant data (Petraitis and Latham 1999, Bertness et al. 2002, Schroder et al. 2005, while theoreticians develop leading indicators of abrupt transitions that demand large amounts of data (Carpenter and Brock 2006, Biggs et al. 2009, Contamin and Ellison 2009. There is little clarity regarding the use of existing data and the design of future studies to detect and mitigate undesired state changes (Bestelmeyer 2006, Groffman et al. 2006. A common, systematic approach to analyzing state changes could allow ecologists to marshal a large body of useful data and detailed knowledge to help society better understand and, ultimately, manage abrupt transitions. Here, we illustrate a general, data-based, and mechanismcentered analysis of abrupt transitions using four datasets from the US Long-Term Ecological Research (LTER) program on pelagic ocean, coastal benthic, polar marine, and semi-arid terrestrial ecosystems. These LTER data include some of the longest time-series available for both causal environmental drivers and biological response variables, and interpretations of associations between the drivers and the response variables are enhanced by experimental and mechanistic studies conducted at the same sites.
We first lay out a synthetic framework for describing abrupt transitions and state changes that can be used to compare and contrast among case studies. We then propose a standard analytical approach that provides strong tests for detecting abrupt transitions between states. This approach revealed unexpected results for the pelagic ocean system for which a ''regime shift'' had been described previously, provided stronger evidence for hypothesized state changes in the coastal benthic ecosystem, and yielded new evidence for state changes in the polar v www.esajournals.org marine and semi-arid terrestrial ecosystems. Our analyses illustrate how to identify and interpret causes of abrupt transitions, and also illustrate limitations common to many datasets used to study abrupt transitions and state changes. We conclude with recommendations for improving ongoing and nascent long-term research programs aimed at detecting and forecasting state changes.

A common framework for describing state change
Studies across a wide range of ecosystems reveal five common data elements used in the recognition and analyses of state change: environmental drivers; triggers; biological responses; response mechanisms; and contextual information ( Fig. 1). We introduce these element categories based on earlier syntheses (Scheffer et al. 2001, Andersen et al. 2009, Suding and Hobbs 2009) and consideration of the datasets presented herein.
State changes in biological responses are caused directly or indirectly by changes in environmental drivers. Drivers are usually abiotic and include changes in climate (e.g., temperature, precipitation), or land-use (e.g., resource extraction, nutrient input rates). Environmental drivers usually are considered ''slow variables'' (e.g., Folke et al. 2004, Carpenter andBrock 2006) because they typically change much more slowly than biological response variables ( Fig. 2A). The textbook example of a slow environmental driver leading to a state change is long-term phosphorus input leading to an abrupt shift from oligotropic (clear blue) to eutrophic (muddy green) lakes (Carpenter and Brock 2006). Drivers can also change abruptly, however, with dramatic effects. Triggers (a.k.a. pulse disturbances) are either abrupt shifts in drivers or singular events, such as droughts, hurricanes, disease outbreaks, invasive species introductions, or fire, that directly affect biological responses (Suding and Hobbs 2009). State changes often are caused by interactions among multiple drivers and triggers (Nystrom et al. 2000, Breshears et al. 2005. Whereas drivers are typically presented as time series concurrent with biological responses ( Fig.  2A), triggers are discrete events in time or relatively short, discrete sections of a time series (e.g., an El Niño period, Holmgren et al. 2006).
Biological responses (a.k.a. response variables or state variables) are used to recognize alternative states (Mantua 2004, Schroder et al. 2005, Andersen et al. 2009). Response variables are especially important because they usually can be measured or monitored easily, and persistent changes in their mean or increases in their variance are used as indicators of state changes. Like environmental drivers, biological response variables typically are represented as time series of the abundance or biomass of individual species or suites of trophically similar species (Daskalov et al. 2007; Fig. 2A).
Response mechanisms describe how drivers and triggers interact and affect responses (Fig. 1). Of particular importance are changes in relationships between drivers and responses caused by positive feedbacks between them that amplify changes in both drivers and responses and reinforce alternative states (Rietkerk et al. 2004). Positive feedbacks often involve complex chains of interactions involving biological and physical v www.esajournals.org processes, including Allee effects (Bourbeau-Lemieux et al. 2011), trophic cascades (Carpenter et al. 1999, Carpenter et al. 2011, habitat fragmentation and extinction cascades (Swift andHannon 2010, He andHubbell 2011), land surface-climate feedbacks (Foley et al. 2003, Cook et al. 2009), or spreading desertification (Peters et al. 2004). Data on response mechanisms are derived most frequently from manipulative experiments, natural history observations, and expert knowledge (Choy et al. 2009).
Finally, contextual information documents characteristics of the environmental setting that can influence driver-response relationships and that can vary among case studies. For example, lake morphometry (Genkai-Kato and Carpenter 2005), stream channel geometry (Heffernan et al. 2008), soil texture , and distance to source populations (Hughes et al. 1999) result in spatial variation in biological responses to drivers and triggers. Similarly, the timing of disturbance events with respect to seasonal period can determine their effects on biological responses (Nystrom et al. 2000). Understanding spatiotemporal context can help to reconcile differences among case studies illustrating general types of transitions and state changes (e.g., Petraitis et al. 2009). Contextual information also can help translate scientific analyses into meaningful policy recommendations and management interventions (Carpenter et al. 2011).

An approach for identifying abrupt transitions and state changes in ecological systems
Three general classes of mechanisms are postulated to produce abrupt transitions: linear tracking, threshold response, and hysteresis Fig. 2. Three classes of driver-response relationships and analytical indicators of transitions and state changes. The top row (A) illustrates time series of driver and response variables in linear tracking, threshold, and hysteretic systems. The second row (B) illustrates that the frequency distribution of the observations should shift from unimodal to bimodal when a threshold or hysteretic change occurs. The third row (C) illustrates how one leading indicator, the variance of the time series, should differ among the three classes of driver-response relationships. As the transition becomes more abrupt and the post-transition state becomes more distinctive from the pre-transition state, the variance should become more peaked at the transition point. The bottom row (D) illustrates changes in the driver-response relationships from linear (in the linear tracking class) to nonlinear (in the threshold class) to hysteretic.
v www.esajournals.org (following Scheffer et al. 2001, Andersen et al. 2009, Suding and Hobbs 2009 (Fig. 2). Note that all three mechanisms can yield patterns that have been referred to as ''thresholds'' in biological response data. An integration of exploratory data analysis, time-series analysis, and linear or nonlinear modeling (see Methods) provide evidence for assigning each case to a class.
The distinction between linear tracking and threshold responses is whether: the distribution of the biological response variable is unimodal vs. weakly bimodal (Fig. 2B); the variance in the biological response is constant vs. increasing slightly as the environment changes or a trigger occurs (Fig. 2C); and the relationship between the environmental driver and the biological response is linear vs. nonlinear (Fig. 2D). Following the terminology of preceding authors, both linear tracking and threshold responses can be reversed; as the driver returns to its initial (pre-change) value, environmental conditions and biological response variables often track them with at most short time-lags. Note that the threshold in ''threshold response'' refers to the nonlinear biological response to a change in driver magnitude, rather than irreversibility.
In contrast to threshold responses, hysteretic responses result from persistent environmental changes, changes in feedbacks between drivers and response variables, or long time lags in biological responses to drivers. In such systems, even if the environmental driver returns to earlier values, the biological response may not return to its earlier state, or does so only slowly, at a markedly different magnitude of the driver, or along a different path from the one it took to reach its new state (Fig. 2D). The functional form of the relationship between environmental driver(s) and biological response(s) typically differs before and after a state change.

Case studies
We examined long-term datasets from four US-LTER programs to characterize abrupt transitions and state changes following our framework; the California Current System (California Current Ecosystem LTER; http://cce.lternet.edu), Southern California Rocky Reef (Santa Barbara Coastal LTER; http://sbc.lternet.edu), Western Antarctic Peninsula (Palmer Station LTER; http://pal. lternet.edu), and Chihuahuan Desert (Jornada Basin LTER; http://jornada-www.nmsu.edu). The California Current System study focused on the abundance of a euphausiid (krill), Nyctiphanes simplex, as a biological response and its relationship to the Pacific Decadal Oscillation Index (PDO) as a primary environmental driver; PDO represents changes in the ocean physical environment that affect krill populations, including advection patterns and water column conditions. The Southern California Rocky Reef study focused on the abundance of the red sea cucumber (Pachythyone rubra) and its relationship to the number of days with large waves (.3.25 m) per year. These large waves disrupt the dominance of P. rubra and allow dominance of macroalgae and associated fauna. The Western Antarctic Peninsula study considered shifts in the abundances of three Pygoscelis penguin species: the Adélie (P. adeliae), chinstrap (P. antarctica), and gentoo (P. papua). These biological responses were considered with respect to changes in the seasonal duration of sea-ice that influences the foraging and breeding biology of these species. Finally, the Chihuahuan Desert study examined changes in the production of the dominant grass species black grama (Bouteloua eriopoda) and its relationship to summer rainfall that governs its production. Details on each case study can be found in the Appendix.

General analytical approach
For each of the four case studies described individually below, we used a sequence of five steps to identify abrupt transitions and characterize state changes with respect to the classes of mechanisms: (1) visualization of temporal patterns in drivers and response variables; (2) locating and statistically testing one or more breakpoints in time-series of response variables; (3) statistical testing of unimodality of frequency distributions of response variables; (4) calculation of temporal variance (a leading indicator used to forecast state transitions) of response variables, and (5) assessment of relationships between response variables and drivers before and after breakpoints identified in (2). Contextual information used to interpret the results was derived from ancillary experimental data, expert knowledge on triggers and response mechanisms, and other v www.esajournals.org natural history information (Appendix). Datasets and associated metadata are archived on, and publicly available from, the Harvard Forest Data Archive (http://harvardforest.fas.harvard.edu/ data/archive.html), dataset HF170. All data manipulation and statistical analyses were performed using base and user-contributed functions in the R language environment (R Development Core Team 2011), as detailed below. The R code used is presented in the Supplement and also available with dataset HF170 described above.
Prior to any analyses, observations of response variables were standardized: z i ¼ (x i Àx )/sd(x). By working in standard deviation units, data and analyses were comparable across the studies. The response and driver variables were unique to each of the four case studies (Figs. 3-6) and time series ranged from 28 to 59 years long. The time series of the responses in each case study included missing data, so modeled values were used in place of missing values. Modeled values were generated from a normal distribution with the mean and variance equal to the running mean and variance, respectively, of the standardized measured values bracketing the missing value(s). For example, in a time series running from 1970-2010, if observations were missing for 1975-1978 and 1980, the modeled values would be sampled from N(mean [z 1974 , z 1979 , z 1981 ], SD[z 1974SD[z , z 1979SD[z , z 1981). Below, we use fz i g to refer to the time series that includes both observed and modeled response variables in standard deviation units.
Temporal patterns in responses were visualized by fitting a locally weighted scatterplot smoother (LOESS) (Cleveland and Devlin 1988) to fz i g. The smoothed curve was fit using the loess function in the R stats library. Default settings were used: a weighted least-squares fit to a fraction of the points in a moving window that spanned three-quarters of the points. The weighting function for each point was proportional to the cube of the distance to each point in the moving window. The curve is fit using a lowdegree polynomial to a subset of the data using a weighted least squares method (Cleveland and Devlin 1988).
Breakpoints in fz i g were identified using the strucchange package (Zeileis et al. 2002). First, the time series was detrended by differencing using the diff function in the R base library. A detrended time series of standardized observations has slope equal to zero, and if there is no breakpoint in the time series, the intercept also would be equal to zero. Breakpoints are years after which the intercept of the detrended time series changes significantly, and detection of one or more breakpoints would suggest that an abrupt transition may have occurred. A combination of three approaches was used to detect breakpoints and to determine the number of breakpoints in the fz i g for each case study. First, a cumulative sum (CUSUM) plot summarized the cumulative sums of differences between each value and the overall mean. A breakpoint was indicated by a sudden change in direction of the CUSUM plot. Because CUSUM plots are ''jagged'' and can indicate many directional changes, residual sums of squares (RSS) and the Bayesian Information Criterion (BIC) were used to identify the number of breakpoints that significantly improved the fit of the CUSUM model (Zeileis et al. 2002). Finally, we examined the statistical significance of each breakpoint identified from RSS and BIC using an F statistic (based on the Chow test statistic, Zeileis et al. 2002). CUSUM, RSS, and BIC plots all indicated either one or two breakpoints in each of the case studies. Because changes in response variables exceeded two standard-deviation units only in the case of the gentoo penguins, however, F statistics were significant only at the a ¼ 0.1 level.
Histograms and density smoothers of fz i g were plotted to determine if the frequency distribution was unimodal or bimodal. Departures from unimodality were tested using Hartigan's dip test (Hartigan and Hartigan 1985) as implemented in the dip function in the R diptest library. This test is very conservative-the distribution of the test statistic is based on asymptotic and empirical samples relative to a uniform distribution. A table of quantiles (P-values) is provided in the file qDiptab in the R diptest library. The power of the test (for a ¼ 0.05) is 80% when sample size ¼ 50; since our sample sizes (excluding missing values) ranged from 27 to 55, we accepted P-values 0.10 as statistically significant evidence for departure from unimodality. The linear tracking model should yield a unimodal distribution of fz i g, whereas a threshold or hysteresis model should yield a bimodal distribution of fz i g (see Fig. 2). (continuation of Fig. 3 legend). Finally, the relationship between the driver and response (E) are illustrated for the initial state (solid symbols, black lines) and post-transition state (open symbols, grey lines). In this figure only, the data from the third state are combined with those from the first state (negative phase of the PDO).  (Carpenter and Brock 2006). We calculated changes in temporal variance of the differenced Fig. 5. Results of analyses for the Western Antarctic Peninsula (A-E), following the same order and rules as for Fig. 3. The three species of penguins are illustrated in three colors (Adélie penguins in black, chinstrap penguins in orange, gentoo penguins in blue). The relationship between the driver and response is shown only for Adélies because there are too few data for the other species.
v www.esajournals.org 9 December 2011 v Volume 2(12) v Article 129 time series using the rollapply function in the R zoo library. The window size used for each case study was the shortest time-interval between breakpoints in the time series; window sizes ranged from seven to 30 years. The temporal variance for years prior to the onset of our moving window could not be calculated (as the number of points available was less than the window size); we indicate those years with dashed lines in Figs. 3D, 4D, 5D, and 6D. We note that using temporal variance as a leading indicator works best for very long time series ()50 observations) of temporally autocorrelated data sampled at high frequency. Real ecological data, such as those analyzed here, are of relatively short duration (,50 observations), and ecologists generally attempt to minimize temporal autocorrelation by sampling less frequently. If the threshold response occurs within the lifespan of the organism, but sampling frequency is on the same time-scale as organism lifespan or generation time (cf. Fig. 7), a state change or threshold response may not be detected. Finally, if observation errors are relatively large or if multiple linear and non-linear processes interact and mute the response variables, changes in variance may not be detected even though state changes have occurred (Scheffer et al. 2009, Brock and Carpenter 2010, Carpenter et al. 2011. Finally, relationships between response and driver variables were examined for the data overall and for data partitioned into before and after breakpoints. For the California Current System data, the data were partitioned into sets when the PDO was either negative (before the first breakpoint and after the second breakpoint) or positive (in between the two breakpoints). For the Southern California Rocky Reef data, we only examined the data before the first breakpoint and after the first, but before the second, breakpoint (only three values for the driver variable were available after the second breakpoint). For the Western Antarctic Peninsula data, we only examined the data for Adélie penguins, because there were too few data for chinstrap or gentoo penguins after their 2004 breakpoints. We used linear (lm) and non-linear (nls) regression in the R stats library to model the relationships between responses and drivers. The expectation for the linear tracking model was that there would be similar response-driver relationships before and after the breakpoint(s), and the expectation for the hysteresis model was that there would be different response-driver relationships before and after the breakpoint(s). For example, a different slope and intercept for a linear regression fitting response-driver relationships or a non-linear versus linear fit for data and after the identified breakpoint would support the hysteresis model Carpenter 2003, Bai et al. 2010).

A pelagic ocean ecosystem: The California Current System
Data collected within the California Current System (CCS) provide an example of abrupt transitions with a linear tracking mechanism (Fig. 3). The CCS includes a major coastal upwelling biome that extends from British Columbia to Baja California. A variety of directional changes in the ocean environment (including rising sea level, oceanic warming, increased density stratification, decreased transparency, acidification, and changes in hypoxia) may be affecting planktonic populations and the pelagic food web. There are also important sources of interannual (e.g., El Niño-Southern Oscillation [ENSO]) and decadal (e.g., Pacific  (Mantua et al. 1997) variability in this ecosystem. Long-term variations in krill abundance are correlated with the PDO (Brinton and Townsend 2003) and time series of N. simplex abundance display abrupt shifts from one persistent state to another, which may imply hysteresis (Fig. 2) and/or a positive feedback mechanism (deYoung et al. 2008). We assessed the evidence for alternative states in the krill population in the southern sector of the California Current where temperate-subarctic, cool-water zooplankton fauna enter from the north, and subtropical, warm-water fauna, including N. simplex (Brinton et al. 1999), enter from the south. This geographic location is therefore likely to be sensitive to changes in large-scale ocean circulation patterns incorporated into the PDO.
The six-decade California Cooperative Oceanic Fisheries Investigations (CalCOFI) record revealed that N. simplex generally was rare when the PDO was in the negative phase (anomalously cool waters in the NE Pacific). Abrupt changes of the PDO from negative to positive were tracked by increases in N. simplex abundance and vice versa (Fig. 3A, B). Strong El Niño (1958-60, 1997 and La Niña (1998-99) events had positive and negative influences, respectively, on N. simplex abundance that interacted with changes in the PDO to accentuate abrupt changes or interrupt relationships with PDO (Appendix).
Although the warm phase between 1977-1998 was a period of consistently high abundance of N. simplex relative to the time period before and after, our data indicate that contrary to previous work (e.g., deYoung et al. 2008) this should not be considered a different ecological ''regime'' and hysteresis is not indicated. The frequency distribution of abundances were unimodal and, most definitively, the relationship between N. simplex abundance (response variable) and the PDO (driver) varied linearly with the variations in the PDO and was identical in both the warm and the cool phases of the PDO. Thus, the California Current System illustrates a case of linear tracking (Hsieh and Ohman 2006), without discrete, definable (or ''preferred'') system states. Such linear tracking may be common in shortlived organisms that can quickly and closely track abrupt changes in drivers.

A coastal benthic ecosystem: The Southern California Rocky Reef
Data from shallow rocky reefs off the coast of Southern California provide evidence of hysteresis due to predation-mediated feedbacks (Fig.  4). The reefs can support either a macroalgaedominated community or one characterized by high densities (.10,000/m 2 ) of the filter-feeding sea cucumber, P. rubra. Spatially extensive sea cucumber-dominated states can persist for decades and dramatically alter reef food webs (Rassweiler et al. 2008, Rassweiler et al. 2010). The loss of macroalgae leads to a reduction in micro-crustaceans and their associated fish predators Schmitt 1989, Schmitt andHolbrook 1990a, b).
Time-series data from nine sites spread along a 5-km stretch of coastline on the north shore of Santa Cruz Island illustrate the mechanisms of abrupt sea cucumber-to-macroalgae transitions ( Fig. 4A; see also Rassweiler et al. 2010). The frequency distribution of annual sea cucumber abundance data revealed evidence of bimodality (Fig. 4C). The first transition from macroalgae to sea cucumber dominance occurred in the late 1980s and was associated with a series of years in which there were few high wave events during winter storms (Fig. 4B). High waves dislodge sea cucumbers from algal beds (Rassweiler et al. 2008), but when winter storms are weak, sea cucumbers competitively displace algae by smothering and killing them.
Frequent, strong storms returned after 1995, but the relationship between days of high waves (driver) and sea cucumber abundance (response) disappeared and sea cucumbers continued to dominate the system (Fig. 4E). Consumption of algal spores by abundant sea cucumbers allowed this species to persist in the face of increased wave disturbance (Rassweiler et al. 2008). This relationship switched to yet another low cucumber state when predatory sea stars colonized the system in late 2002 (Appendix). Thus, this case conforms to a hysteresis model in which stabilizing feedbacks conferred resilience with respect to the environmental driver.
A polar marine ecosystem: The Western Antarctic Peninsula The Western Antarctic Peninsula (WAP) provides another example of hysteresis due to the v www.esajournals.org effects of multiple, interacting drivers (Fig. 5). Since 1950, annual mean air temperature in some regions has increased by 28C, and winter air temperature has increased by nearly 68C (Smith et al. 1996, Vaughan et al. 2003) These climatic changes have caused longterm reductions in the regional extent and duration of winter sea-ice (Smith and Stammerjohn 2001, Stammerjohn et al. 2008a, Stammerjohn et al. 2008b), a proximate driver of directional environmental change in the WAP marine ecosystem (Smith et al. 2003, Ducklow et al. 2007, Moline et al. 2008. One important change is a poleward shift in breeding ranges of three closely related penguin species; the Adélie, chinstrap, and gentoo (Forcada et al. 2006, Ducklow et al. 2007, Forcada and Trathan 2009). There is considerable debate regarding the environmental drivers of change in penguin breeding population dynamics (Patterson et al. 2003, Forcada and Trathan 2009, Trivelpiece et al. 2011. Nearly 40 years of data collected from the Palmer Archipelago near Anvers Island, Antarctica (Appendix: Fig. A2) illustrate abrupt declines in the Adélie penguin breeding population beginning in 1993, and abrupt increases in numbers of breeding chinstrap and gentoo penguins beginning in 2004 (Fig. 5A). Bimodality of annual abundance data was not clearly evident (Fig. 5C), nor did temporal variance (Fig. 5D) illustrate dramatic changes before or during the observed population changes. However, analysis of the relationship between the proximate driver (sea-ice duration; Fig. 5B) and Adélie penguin breeding population size revealed that prior to the 1993 breakpoint, the abundance of this species was essentially unresponsive to variation in sea-ice duration, however after 1993 these variables were strongly and positively correlated (Fig. 5E), conforming to the hysteresis model. We did not examine driverresponse relationships for chinstrap or gentoo penguins because only five data points on yearly numbers of breeding pairs have been obtained since the 2004 breakpoint. Progressive climate warming resulted in an abrupt transition operating through multiple, cascading ecological drivers and feedbacks, including reduced sea-ice duration, changes in terrestrial snowfall accumulation that affect penguin breeding biology, and feedbacks between Adélie population reductions and predator efficiency (Appendix).

A semi-arid grassland ecosystem: The Chihuahuan Desert
Data from northern Chihuahuan Desert grasslands provide an example of hysteresis involving a strong trigger and novel feedbacks (Fig. 6). These grasslands were dominated historically by black grama grass (Bouteloua eriopoda), but during the last 150 years, black grama grasslands have shifted to shrublands dominated by xerophytic woody plants. Similar shifts from grasslands to shrublands have occurred in semi-arid systems throughout the world (Archer 1995). Historically, black grama grass persisted through episodic droughts, and shrub cover within black grama grasslands was limited by competition for water, limited shrub seed dispersal, and possibly periodic fire (Peters and Gibbens 2006). Heavy cattle grazing on black grama grass during drought periods is believed to have initiated the grassland-to-shrubland transition (Appendix). It has not been clear, however, how rapidly the initial grassland loss takes place and therefore how best to employ monitoring strategies to prevent it (Bestelmeyer 2006).
Time-series data on annual production of black grama grass collected during the mid-1900s from two pastures in the Jornada Experimental Range near Las Cruces, NM, USA, indicate the start of an abrupt transition in 1948. In that year, there was no black grama production (Fig. 6A), and this lack of production coincided with the onset of a prolonged drought (Fig. 6B). Several lines of evidence suggest that this system conforms to the hysteresis model. First, annual production was bimodal (Fig. 6C), indicating two alternative states. Second, black grama production exhibited an increase in temporal variance during the transition (Fig. 6D) associated with a period of low and variable growing-season (July-September) rainfall (Fig. 6B). Third, driver-response regressions show that prior to 1948, black grama production had a positive relationship to growing-season precipitation (Nelson 1934). After 1948, however, this relationship weakened and overall production was low regardless of growing-season rainfall (Fig. 6E). The shift in black grama production was very abrupt, never attaining previous high values after 1950. A positive feedback between soil erosion and low grass cover appears to have precluded grassland recovery after a return to higher levels of precipitation (Appendix).
To summarize, we used a common approach to determine if and when abrupt transitions occurred, to evaluate leading indicators that could forecast the transitions, and to match each case to the appropriate class of mechanisms (Figs. 3-6). The timing of abrupt transitions was successfully identified in all four cases. An increase in variance that could serve as a leading indicator was observed only in the Chihuahuan Desert case due to the extreme interannual fluctuation preceding grassland collapse. The linear tracking model was indicated for the California Current case due to unimodality in the distribution of biological response values and linearity in the driver-response relationship. The hysteresis model was indicated in the other three cases due to varying combinations of evidence, including the strong nonlinearity in the driverresponse relationship for the West Antarctic Peninsula, and both bimodality of biological responses and nonlinear driver-response relationships in the Southern California Rocky Reef and Chihuahuan Desert cases. In all four cases, plausible response mechanisms supported the classification of the case to the general mechanism.

A common approach
These case studies illustrate that abrupt transitions and state changes not only can be identified, but also can be understood via a suite of general concepts (Fig. 1) and relatively simple methods. Although the availability of longrunning time series of both drivers and responses has been limited (Andersen et al. 2009, Carpenter et al. 2011, long-term data now can be accessed from LTER and related sites (http://ecotrends. info), and many institutions worldwide are investing considerable resources establishing new ecological observation networks (e.g., National Ecological Observatory Network, Global Lake Ecological Observatory Network, Ocean Observatory Initiative, Paleoecological Observatory Network). The sequence of methods used here, including an objective evaluation of abrupt transitions and alternative states, examination of frequency distributions of response variables, consideration of variance patterns used in forecasting, and analysis of patterns and mechanisms of driver-response relationships can be applied to many of these datasets.
A common, systematic approach applied across different datasets will advance a general understanding of abrupt transitions and state changes. Such a common approach is especially important now, as abrupt, often irreversible transitions are forecast to increase as climatic change accelerates (Millennium Ecosystem Assessment 2005, Fagre et al. 2009); a coherent, integrated strategy is needed to manage and mitigate the expected state changes. Our comparative exploration of case studies suggests some lessons for future analyses of existing data and guidance for new observation and monitoring networks embarking on long-term studies.

Leading indicators may have limited utility
Mathematical modeling (Scheffer et al. 2009) and empirical studies conducted in temperate lake ecosystems (Carpenter et al. 2011) suggest that increased variance in the time series of biological responses can be used to forecast abrupt transitions. Many systems, however, may show no change or even decreases in variance (Hastings and Wysham 2010). Our analysis of four different systems, three of which showed clear hysteretic patterns between drivers and responses, suggest that this leading indicator must be carefully scaled to the time-scale of dynamics in the biological response variable (e.g., organism lifespan; Fig. 7). Short-lived organisms can track abrupt changes in drivers closely. Thus, some transitions, such as those in the California Current krill (Fig. 3), may appear abrupt until rescaled to the short lifespan of this organism (Fig. 7). In contrast, especially when lifespan matches the dominant time scale of environmental variability (Hsieh and Ohman 2006), other transitions may appear gradual, but actually occur quite abruptly relative to the organism's lifespan (e.g., penguins: Fig. 5A; black grama grass: Fig. 6A).
To be informative, leading indicators of rising variance require many highly autocorrelated samples collected within the lifespan of the sentinel organism of interest. The traditional v www.esajournals.org ecological emphasis on temporal trend as opposed to variance has led researchers to undervalue the importance of regularly sampled time series and fine temporal intervals. Detection of abrupt transitions and state changes require time series without missing values (or ones that can be reasonably modeled). If sampling designs capture weakly-or un-correlated measures of abundance, or if studied organisms are longlived and transitions occur rapidly (i.e., between samples), measures of temporal variance may not be informative. In such cases, it would be worthwhile to identify faster-changing variables (Carpenter et al. 2011) that reflect organismal performance in populations, such as physiological status. However, if transitions are not caused by ''slow'' variables but instead are caused by abrupt, unpredictable triggers acting within vulnerable systems, variance-based leading indicators may provide only limited information (Hastings and Wysham 2010). In such cases, mean values of a slow variable might be used to signal heightened susceptibility to abrupt transition given a trigger, but the transition itself may not be predictable.

Driver-response relationships are powerful tools
Researchers should hesitate to infer response mechanisms based solely on the presence of threshold patterns in biological response variables; analysis of driver-response relationships provide stronger tests of such inferences. For example, the linear tracking model (Fig. 2) may appear to have abrupt transitions when biological responses track abrupt changes in drivers, as in the California Current System (Fig. 3A). Such observations have been used to suggest the existence of alternative stable states (deYoung et al. 2008) that is not supported by our analysis (Fig. 3E). Driver-response relationships can then be explained with regard to specific response mechanisms (e.g., grass cover loss leading to soil erosion feedbacks in the Chihuahuan Desert case; Appendix).

Context is critical
In the four cases that we examined, the historical context and the location in which the study was conducted relative to physical processes occurring at larger spatial scales both had important consequences for the observed dy-namics. The location of the California Current System and Western Antarctic Peninsula studies relative to biogeographic transition zones, the Chihuahuan Desert study on sandy soils during a period of comparatively high grazing pressure, and the Southern California Rocky Reef study area relative to the shifting southern range limit of sea-stars each influenced the patterns observed in their respective time-series of biological responses (Appendix). Data on the same response variables collected at other locations or times might yield different results or reveal how large-scale forces such as ocean circulation, regional climate, physiography, or soils mediate abrupt transitions (Rietkerk et al. 2004, Williams et al. 2011. As case studies of abrupt transitions accumulate, researchers should ensure that spatial and temporal measurement scales of drivers and response variables are recorded and are appropriately congruent. Researchers should also be alert for changes in context when comparing studies.

Multiple lines of evidence aid interpretation
Different kinds of state changes were identified by different analyses. Interpretations of state changes were greatly facilitated by consideration of data and observations apart from the driver and response time series. The choice of the ''right'' driver and the ''right'' response variable was based on detailed short-and long-term experiments along with ancillary information and anecdotal observations that provided important clues to the interpretation of time-series data. In spite of the case-specific nature of response mechanisms, we predict that a systematic review of additional cases will reveal a limited set of classes of interactions between drivers, triggers, and responses (cf. Fig. 2). This framework can guide future investigations and promote a datasupported understanding of abrupt transitions (Walker and Meyers 2004).

Can we manage state changes?
Long-term, multi-faceted case studies and datasets can provide retrospective explanations of transitions and state changes for specific cases, but can they provide useful information for proactive management? We suggest that a body of transition datasets representing different kinds of ecosystems would provide several useful v www.esajournals.org insights for managing state changes. First, they would suggest the appropriate drivers, triggers, and biological response variables to be emphasized in monitoring, as well as the spatial and temporal design elements needed for detection. Second, depending upon the drivers, such analyses would indicate useful strategies for managing the effects of drivers. In some cases, such as when response variables are affected by broad-scale physical drivers (e.g., California Current System) direct control is not feasible but forecasting and adaptation approaches could be developed. Third, case studies can be used to evaluate the abruptness of transitions, particularly with respect to organisms' lifespans (Fig. 7). Cases of high abruptness, such as in the semiarid grasslands and Antarctic peninsula, suggest that management reactions to indicators of transition must be especially rapid (e.g., adjust stocking rates in drought periods and establish institutions that can respond rapidly) (Meze-Hausken et al. 2009). Finally, case studies can be used to indicate the potential utility of early warning indicators to forecast transitions given their abruptness and the feasibility of measuring the appropriate attributes with sufficient temporal and/or spatial resolution. Many abrupt transitions may ultimately need to be managed according to a precautionary principle that acknowledges our limited ability to develop indicators of imminent transition or to respond rapidly enough to such indicators (Contamin and Ellison 2009). In those situations, case studies might indicate simple values in driver or biological response levels that are related to the likelihood of abrupt transition. To make informed choices among the wide range of possible strategies for detecting and managing abrupt transitions, ecologists and policymakers must commit to sustaining, renewing, and initiating observational platforms in multiple sentinel sites. The resulting data can, as we have shown, produce useful maps for navigating our changing world.

Detail on the California Current System case study
To assess the likely mechanism underlying the relationship between variations in the California Current euphausiid (krill) Nyctiphanes simplex and the Pacific Decadal Oscillation (PDO) first shown by Brinton and Townsend (2003) and reexamined in this paper (Figs. 3A-E), we analyzed the life-history structure of N. simplex from 1951 to 2009. Zooplankton were sampled in the upper 210 m or 140 m of the water column using 0.5-mm mesh plankton nets (Ohman and Smith 1995) and analyzed by E. Brinton. The stations used in the analysis were from Southern California (Lavaniegos and Ohman 2007); the California Cooperative Oceanic Fisheries Investigations (CalCOFI) lines 80 through 93, from shore to station 70, springtime cruises, and nighttime samples only. This station pattern differs slightly from that used by Brinton and Townsend (2003). Annual averages of the Pacific Decadal Oscillation (Mantua et al. 1997) were obtained from the monthly values posted at: http://jisao. washington.edu/pdo/PDO.latest.
Counts of krill during these 48 years were available for four life-history stages: calyptopis, furcilia, juveniles, and adults. The proportional composition of each stage (Fig. A1) was used to differentiate between two primary mechanisms by which changes in the physical environment, as represented by the PDO, might have influenced euphausiid abundance: altered advection of organisms into or out of the study region, and altered in situ changes in water column conditions (e.g., temperature, food, predators) that can affect population growth of N. simplex.
The primary breeding center of N. simplex is located off Baja California (Brinton et al. 1999), to the south of our study site. Abrupt increases in advection of organisms into our region from the south, reflecting individuals introduced from a population showing regular recruitment in a favorable habitat, would have rapidly increased the proportion of younger individuals (calyptopis larvae) in the population. Conversely, if conditions changed favorably in situ, without corresponding changes in advection, the contribution of calyptopis larvae would have increased more gradually over time. Finally, if in situ conditions for N. simplex deteriorated without a corresponding change in advection, there should have been a decline in relative abundance of larval stages due to reduced egg production by adults.
Inspection of Fig. A1 relative to the two abrupt transitions identified in the population time series (Fig. 3A) suggested that altered advection was the more plausible hypothesis. The 1976-77 abrupt increase in total abundance of N. simplex was accompanied by a simultaneous increase in the proportion of larvae, without a temporally lagged response. The explanatory power of the PDO for the temporal variability in N. simplex is corroborated by modeling variations in N. simplex as an autoregressive (AR-1) process related only to present and one previous state of the PDO, which shows excellent agreement with observations. Furthermore, the rapid decline in total abundance in 1999-2000 was not accompanied by a gradual diminution in contribution by larvae. Rather, the proportion of larvae remained roughly constant, although production was intermittent. Hence, we conclude that the predominant mechanism underlying rapid changes in the euphausiid-PDO time series was the introduction or flushing out of individuals through altered transport. Once introduced into the study site, N. simplex were able to reproduce and survive for extended periods of time because of more favorable conditions in situ, but the rapid increases/decreases in abundance were initiated by altered advection. Further support for this interpretation comes from observations of responses by N. simplex to major El Niño events. These events typically have resulted only in transient increases in abundance of N. simplex (Fig. 3A, B; see also Brinton and Townsend 2003) initiated by transport from the south. The unusual El Niño of 2009-2010 was not accompanied by changes in N. simplex abundance because this particular event propagated through atmospheric teleconnections rather than through altered ocean advection (Todd et al. 2011).

Detail on the Southern California Rocky Reef case study
We used data on abundance of red sea cucumbers (Pachythyone rubra) from nine sites on the north shore of Santa Cruz Island, CA. Sites were situated along a 5 km stretch of coastline (centered on 34.058N, 119.7378W), with six sites established in 1982 and three more added in 1989 (for a complete description of the sites see Holbrook 1986, Holbrook andSchmitt 1989). Sampling was usually annual between September and November; in some instances not all sites were sampled every year. Sites were similar in depth, slope, exposure and initial benthic community structure. At each site there were two fixed 40 m transects, one each along the 6-m and 9-m isobath. The percent cover of P. rubra, understory macroalgae (mainly species in the genera Eisenia, Laurencia, Gelidium, Rhodymenia, Codium and Corallina), sea urchins, and Macrocystis pyrifera were estimated using point-contact methods (eight randomly located points per meter per transect). The sea-star Pycnopodia helianthoides was counted in a 2-m wide swath centered on each transect. Abundance of P. rubra was estimated as the average percent cover across all transects sampled in each year. Similarly, abundance of P. helianthoides was estimated as the mean density across all transects sampled in each year.
Data on wave heights were taken from buoys operated by the National Oceanic and Atmospheric Administration (www.ndbc.noaa.gov). Because no single buoy operated without interruption between 1982 and 2008, data from three buoys were combined. The East Santa Barbara buoy (No. 46053) is located nearest to the sites where organismal cover and abundance data were collected (23 km NNW of the study sites) and, therefore, wave heights recorded at this buoy were used when available. On days when data were not available from the East Santa Barbara buoy, data from the Point Arguello and Santa Maria buoys (Nos. 46023 and 46011, ;135 km NW of the monitoring sites) were used to estimate wave height and water temperature in the east channel. Estimates were made based on linear regressions developed from days when all three buoys were operational (Rassweiler et al. 2008).
To quantify wave intensity, we calculated the number of days each year when maximum significant wave height exceeded 3.25 m. October 1 was the cut-off between years, because biological sampling typically occurred near this date. We chose 3.25 m as our definition of a large storm because previous experiments revealed that the competitive effects of algae on P. rubra abundance only occurred when waves exceeded this height (Rassweiler et al. 2008, Rassweiler et al. 2010). We did not include wave data from the summer period of each year (May through September), because summer swells typically come from a southerly direction and the northern shores of the Channel Islands, where our sites are located, are sheltered from these waves.
One challenge in analyzing state changes is that there can be more than two states. In this case study, an exclusive focus on P. rubra initially suggested only two states: prior to 1987 the sea cucumbers were nearly absent, from 1987 to 2002 they were very abundant although variable, and after 2002 they returned to their low density state, approaching the densities initially observed in the early 1980s (Fig. 4A, E). However, when other invertebrates were considered, it became apparent that the post-2002 low abundance state was not equivalent to the pre-1987 state. Rather, the post-2002 system consists of different species and is maintained by a new mechanismpredation. The third state was initiated by the sudden arrival of the predatory sea star, P. helianthoides, a mobile and voracious predator, with well-established potential to decimate echinoderm populations (Mauzey et al. 1968, Pearse andHines 1987). In 2003 when sea stars first became abundant they were typically large-35 cm or more across-suggesting that they immigrated into these sites, either from deeper water or from the western end of the island where they have been previously observed (Eckert 2007).
The role of sea stars in the initiation and maintenance of the third state illustrates that an interaction between multiple environmental factors triggered the post-1987 shift into the high P. rubra state. The absence of predators alone was not sufficient to have caused this shift; before 1987, P. rubra was rare even though predators were absent. Similarly, it is unlikely that low waves could have triggered a shift into the high density state if P. helianthoides had been present, because the sea stars exert strong top-down control on P. rubra abundance. Our results underscore the complex nature of state changes. Different mechanisms can initiate, maintain, or end a state, and interactions between multiple drivers may be necessary to trigger shifts in states.

Detail on the Western Antarctic Peninsula case study
Species comprising polar marine systems have evolved life histories associated with the presence or absence of sea-ice, often broadly termed sea-ice obligate or sea-ice intolerant species, respectively (Ducklow et al. 2007, Moline et al. 2008, Siniff et al. 2008. Pygoscelis penguins of the Western Antarctic Peninsula (WAP, Fig. A2) integrate environmental variability over large spatio-temporal scales due to their longevity and spatially extensive foraging (Fraser and Trivelpiece 1996). Relationships between environmental drivers and penguin population dynamics (Fig. 5A, B, E) reflect life history integration of this variability, and the abundance and distribution of these species provided some of the earliest evidence of rapid climate-induced change in the WAP (Fraser et al. 1992, Woehler et al. 2001, Forcada et al. 2006, Siniff et al. 2008, Forcada and Trathan 2009, Gorman et al. 2010, Trathan et al. 2011, Trivelpiece et al. 2011. Physical oceanographic processes occurring along the WAP are important proximate drivers of changes in regional climatology (Thompson and Solomon 2002, Marshall et al. 2004, Ducklow et al. 2007, Martinson et al. 2008, Meredith et al. 2008. Interactions between climate phases and physical oceanography has resulted in displacement of the cold, dry polar climate that historically dominated the region by a warm, moist sub-Antarctic system characteristic of the northern WAP and Scotia Arc (Ducklow et al. 2007).
Penguin population data in this case study span nearly four decades, a period during which sea-ice extent decreased by 50% and sea-ice duration decreased by 85 days (Smith et al. 2003, Stammerjohn et al. 2008b). Number of breeding pairs of Pygoscelis penguins has been estimated annually since the mid-1970s from surveys of nesting individuals on seven islands within 15 km of Palmer Station, a US scientific research station located on Anvers Island (Fig.  A2). Most of the data used in these analyses were based on numbers obtained immediately following peak clutch completion (November-December). In the few years where this peak was missed due to weather and sea-ice conditions hindering island access, the next survey conducted closest to this period was used. During 1980During , 1984During , 1985 and 1988, regional totals were estimated from partial surveys (i.e., data from islands not surveyed were estimated based on percent increases or decreases on adjacent islands that were surveyed). Analyses to examine relationships between sea-ice duration and Adélie penguin population response were lagged by four years to account for delayed reproductive maturity of these species (Ainley 2002), however, results were qualitatively similar for lags equal to zero and five. Following (Stammerjohn et al. 2008b), annual sea-ice duration was based on the number of days that elapsed between the first day of advance and the first day of retreat for the Palmer LTER study region near Anvers Island (Fig. A2); an ''ice year'' begins in mid-February of year y and ended in mid-February of year y þ 1.
Since 1975, the breeding population of the true Antarctic, sea-ice obligate, Adélie penguin (Ainley 2002) along the Palmer Archipelago has declined by 85% (Fig. 5A). The breakpoint in Adélie population dynamics occurred in 1993 (Fig. 5A); this response is temporally consistent, given the species lag in reproductive maturity, with the poorest sea-ice conditions evident in the remote sensing record Hofmann 2003, Stammerjohn et al. 2008b) and the lowest abundance of Antarctic krill (Euphausia superba) from tables in these reports for pastures 2 and 9, which were dominated by this species. Estimates were based on an annually varying number of 15-m long 3 10-cm wide transects. Transects were added until the standard error of the estimate was within 10% of mean production value. On each transect, 100 plants were measured and the height of grazed and ungrazed tillers was recorded. Standing crop of different perennial grass species was estimated by clipping all aboveground grass parts, air-drying them, and weighing them (Jornada Forage Crop Report, 1942, Jornada archives). Areas of each pasture that were not dominated by black grama grass (due to variation in soils) were excluded from sampling. A utilization scale (Lommasson and Jensen 1943) was used to estimate the percent of grazing use for each species, which was averaged over hundreds of plants (Fig. A3). Utilization values equaled the percentage of the recommended biomass removed (35% at that time), and were determined from a ''large number of transects'' randomly placed throughout each pasture in each year; values over 100% (i.e., more than 100% of the 35% recommended use) indicated overgrazing. Precipitation data were from the West Well rain gauge of the USDA Rain Gauge Network (http://usda-ars.nmsu.edu/ data_long-term-datasets.html), which lay at the southwest and northwest corners of pastures 2 and 9, respectively. We paired growth year black grama grass production (measured in fall) with the monsoonal rainfall totals (July-September) of that same year. Known limitations of the data include: (1) a lack of precise spatial relationships between production, patchy rainfall, and pasture utilization; and (2) complex relationships between intra-and interannual rainfall and plant production that are not reflected in the data (Schwinning andSala 2004, Yahdjian andSala 2006).
The state change data reflect the loss of black grama grass cover in discrete areas (i.e., stage I of the ''multi-stage model''). The observed state change was initiated by severe drought years occurring during a period with intermittently high levels of utilization (cattle grazing; Fig. A3). Years of very low summer (July-September) rainfall led to years of relatively low black grama production (Fig. 6A, B, E). The failure to reduce livestock numbers during those years led to overgrazing, and successive years of overgrazing (measured as high utilization values) were followed by years of reduced black grama grass production (Fig. A3). That dry, extremely windy conditions occurring in low grass cover conditions could initiate extensive, severe soil erosion and subsequent collapse of black grama grassland became widely appreciated in the early 1950s (see also Okin et al. 2006). These data, however, reveal how rapidly (over 2 years) these effects can lead to persistent reductions in black grama.