Ecological thresholds are characterized in the literature as abrupt changes in a response variable as a consequence of continuous change in an independent variable (Muradian 2001, Groffman et al. 2006, Dodds et al. 2010). The threshold is defined as the point where this function shows a change in the value (or sign) of the slope. Theoretical and empirical support exists for the hypothesis that some ecosystems show abrupt, nonlinear responses to natural or anthropogenic disturbance (May 1977, Connell and Sousa 1983, Bellwood et al. 2004). However, estimating the location of an ecological threshold or predicting when an ecosystem is approaching a threshold is challenging (Scheffer et al. 2009). A variety of statistical and mathematical approaches have been developed to detect ecological thresholds, and development of these approaches is an active area of research (Dodds et al. 2010, King and Baker 2010).

The traditional application of threshold models in aquatic ecology has been to characterize stressor–response relationships, such as the responses of primary producers to nutrients (Qian et al. 2003), threatened species to habitat loss (Huggett 2005, Hilderbrand et al. 2010), or coral reefs to global change (Bellwood et al. 2004). However, recovery trajectories following removal of a stressor also can be nonlinear (Lake et al. 2007), in part because other sources of disturbance or regional climatic variation might continue to affect recovery. Community responses to multiple perturbations are rarely additive (Paine et al. 1998), so superimposing new disturbances on a chronically stressed system might significantly delay recovery and could cause a hysteresis response (Scheffer et al. 2001, Bellwood et al. 2004). In situations where recovery trajectories are nonlinear or hysteretic (Groffman et al. 2006), we suggest that threshold models can be used to characterize responses after the removal of a stressor. Just as thresholds of resistance are estimated by the relationship between ecological responses and stressor levels, recovery thresholds are determined by the relationship between ecological responses and time since stressor removal.

A better understanding of recovery thresholds and nonlinear recovery trajectories could improve our ability to assess restoration effectiveness in aquatic ecosystems and ensure that restoration efforts are targeted for those systems where improvements are most likely (Poole et al. 2004). A recent survey conducted by the US Environmental Protection Agency (EPA) estimated that ∼⅓ of the streams and rivers in the contiguous US are significantly degraded or polluted (USEPA 2000). Expenditures for restoration of these systems, which averaged ∼1 billion USD/y from 1990 to 2003, are expected to increase, and most of these funds will be used for improving water quality (Bernhardt et al. 2005). Despite these significant investments, our ability to assess effectiveness of restoration of contaminated ecosystems is seriously limited by inadequate study designs (Bernhardt et al. 2005), failure to consider ecological theory (Lake et al. 2007), and a lack of agreement regarding appropriate measures of recovery (National Research Council 2007). Developing objective criteria that define restoration success is challenging because apparent recovery of some indicators, such as abundance or species richness, can occur despite persistent alteration in community composition and loss of ecological resilience (Berumen and Pratchett 2006, King and Baker 2010). The Humpty–Dumpty model of recovery (Pimm 1991) predicts that some disturbed ecosystems might never return to predisturbance conditions, or that recovery might not occur until stressor levels are reduced far below those that triggered the initial response (Scheffer et al. 2001). An improved understanding of these complex recovery trajectories might allow water resource managers to estimate the level of stressors that must be achieved to restore aquatic ecosystems following perturbation.

High elevation streams in the southern Rocky Mountain ecoregion provide a unique opportunity to investigate long-term recovery from water-quality degradation in systems simultaneously experiencing effects of global change. Heavy metal contamination from historical mining operations is a major environmental stressor in many of these streams (Clements et al. 2000, Niyogi et al. 2001). A survey of 79 randomly selected streams in the Southern Rocky Mountain Ecoregion showed that ∼23% of these sites were degraded by heavy metals (Clements et al. 2000). Superimposed on the effects of metals are hydrological alterations (Barnett et al. 2005, Mote et al. 2005), increased atmospheric deposition of NO_x and SO₄ (Baron 2000), and increased exposure to ultraviolet radiation (Clements et al. 2008). Interactions among these stressors over the next several decades probably will complicate our ability to assess recovery from heavy metal contamination.

We used results of a long-term research program on the upper Arkansas River, a metal-polluted stream in Colorado, to investigate threshold responses of benthic macroinvertebrate communities to improvements in water quality. Over the past 17 y (1989–2006), we have measured physicochemical characteristics, habitat quality, heavy metal concentrations, and macroinvertebrate community structure seasonally (spring and autumn) at several locations upstream and downstream from California Gulch, a US EPA Superfund Site. In 1992, state and federal natural resource agencies initiated a large-scale restoration program designed to improve water quality and reestablish trout populations in the watershed. Major restoration activities in California Gulch were conducted between 1992 and 1999 and included construction of water treatment facilities (1992–1994) and removal of mine tailings (1994–1999). We applied a new threshold approach, the significant zero crossings model (SiZer), to evaluate long-term changes in benthic macroinvertebrate communities after completion of restoration and to determine whether recovery thresholds exist.

Methods

SiZer is a relatively new statistical approach that characterizes multiple ecological thresholds based on significant changes in the slope of a relationship between state variables and predictor variables. Details of the SiZer approach and its application to ecological thresholds have been published previously (Chaudhuri and Marron 2000, Sonderegger et al. 2009). Briefly, the SiZer model applies a nonparametric smoother to a predictor–response relationship and then examines the derivatives of the smoothed curve to identify a threshold. SiZer then defines an ecological threshold as a point where the 1^st (or n^th) derivative changes significantly. For example, using a nonparametric smoother, such as locally weighted polynomial regression (LWPR; Fan and Gijbels 1996), we can obtain an estimated response function, its derivatives, and 95% confidence intervals (CIs). A tuning parameter (the bandwidth [h] in LWPR) controls the smoothness of the curve, which theoretically can vary from linear regression (which assumes no threshold) to fitting the data perfectly. All points along the x-axis are then classified into 1 of 3 states based on whether the estimated 1^st derivative of the relationship is: 1) significantly positive (i.e., CI of the 1^st derivative contains only positive values), 2) significantly negative (the CI contains only negative values), or 3) possibly 0 (the CI contains 0). Thresholds are defined as significant (p < 0.05) changes in the state of the derivative and are displayed as a function of h.

A unique aspect of the SiZer approach is that it examines changes in the derivative as a function of a predictor variable across all reasonable values of h and displays this information in a single image, the SiZer map. An investigator can use the SiZer map to evaluate support for a threshold relationship across a variety of smoothing levels. The different values for h essentially reflect patterns in the data viewed at different levels of resolution of the predictor variable. Another unique aspect of SiZer is that it easily handles multiple thresholds and displays graphically all significant changes along the predictor variable.

To illustrate the application of SiZer to identify ecological thresholds, we present 2 analyses using long-term water-quality and macroinvertebrate data collected from Arkansas River station AR3 (Clements 2004), ∼100 m downstream from the California Gulch Superfund Site. We collected quantitative macroinvertebrate samples (n  =  5) seasonally (spring and autumn) with a 0.1-m² Hess sampler (350-µm mesh net) from shallow riffle areas. We identified most organisms (except chironomids) to either genus or species and chironomids to subfamily. We collected water samples for metals analysis on each sampling occasion and determined concentrations of Zn, the predominant metal in the system, with atomic absorption spectrophotometry.

First, we used canonical discriminant analysis to examine separation and overlap of macroinvertebrate communities among years based on abundance of the 20 dominant taxa (1989–2006) at this station. We described long-term changes in abundance of metal-sensitive taxa, including total number of heptageniid mayflies, Paraleptophlebia sp. (Ephemeroptera∶Leptophlebiidae), Heterlimnius corpulentus (Coleoptera∶Elmidae), and Pericoma sp. (Diptera∶Psychodidae). We identified these metal-sensitive taxa based on results of long-term monitoring in the Arkansas River and mesocosm experiments conducted with benthic communities (Clements 2004). Second, we used SiZer to identify ecological thresholds in multivariate space based on long-term changes in community composition by examining the relationship between the 1^st canonical axis and year.

Results

Concentrations of Zn declined after completion of restoration activities in California Gulch, but this trend was complicated by significant seasonal and annual variation associated with stream discharge (Fig. 1). Zn concentrations generally were elevated during spring runoff, especially in years with high stream discharge (e.g., 1995–1997). Consistent with these long-term improvements in water quality, abundance of metal-sensitive taxa increased from 1989–2006 and approached the reference station value in 2002 (Fig. 2A). At h  =  1.5, which is a reasonable representation of these data, the slope of the relationship between abundance and year changed in several places. The associated SiZer map for these data (Fig. 2B) indicates how the 1^st derivative of this relationship changed significantly at different values of h (note that h is expressed on a log₁₀ scale). Reading from left to right at a specific value of h, the shading changes indicate where the slope parameter transitions from “significantly different from 0” to “not significantly different from 0” (or vice versa) at a 95% CI. At h  =  1.5 (log₁₀h  =  0.18), the 1^st derivative changed from significantly positive to 0 in 1993, from 0 back to significantly positive in ∼1996, remained positive until ∼2002, and then changed back to 0. We define these significant transitions as ecological thresholds.

Figure 1

Long-term (1989–2006) variation in stream discharge (bars) and Zn concentration (line) at Arkansas River station AR3, ∼100 m downstream from the California Gulch Superfund Site. The 2 arrows indicate initiation and completion of major restoration activities at the California Gulch Superfund Site. Seasonal variation in Zn concentration is illustrated by showing autumn samples as solid circles and spring samples as open circles.

Figure 2

A.—Locally weighted polynomial regression analyses (with 95% confidence intervals) showing the relationship between abundance of metal-sensitive taxa and year at station AR3. The dashed horizontal line shows mean abundance of metal-sensitive taxa at station AR1, the upstream reference site, calculated using macroinvertebrate data collected from 2000 to 2006. B.—SiZer map showing changes in the 1^st derivative as a function of year and bandwidth (h). A significantly increasing 1^st derivative is shown in light grey and a possible 0 derivative is shown in dark grey. The horizontal distance between the 2 white curves ( =  2h) represents the effective width of the locally weighted polynomial. The horizontal dashed line corresponds to h  =  1.5 (log₁₀h  =  0.18), and the arrows indicate significant thresholds at this value of h.

Results of multivariate analysis based on abundance of the 20 dominant taxa showed alternating patterns of gradual change followed by sudden transitions in community composition between 1989 and 2006 (Fig. 3A). Results of LWPR analysis (h  =  1.5) for the relationship between canonical variable 1 (which explained 58% of the variation) and year revealed a highly complex function with several instances where benthic communities changed over short time periods (Fig. 3B). The corresponding SiZer map for these data indicated several significant changes in the 1^st derivative between 1993 and 1997 and reflected abrupt shifts in community composition during this period (Fig. 3C). Benthic communities changed steadily from 1999 to 2002, and a final threshold (from significantly positive to 0) was observed in 2004, after which community composition remained stable for the duration of the study.

Figure 3

A.—Mean (±1 SE) values for canonical discriminant variables 1 and 2 by year (1989–2006) and % of variation explained by each variable in a canonical discriminant analysis of the Arkansas River data based on abundance of the 20 dominant taxa. B.—Locally weighted polynomial regression analysis (bandwidth [h]  =  1.5) showing the relationship between canonical variable 1 and year. C.—SiZer map showing changes in the 1^st derivative as a function of year and h. The color scheme is similar to Fig. 2 except that a significantly decreasing 1^st derivative is shown in black. The horizontal distance between the 2 white curves ( =  2h) shows the effective width of the locally weighted polynomial.

Discussion

Our results suggest that the ecological threshold concept can be extended beyond conventional stressor–response relationships to provide important insights into recovery processes following the removal of a stressor. In the Arkansas River, total abundance of metal-sensitive species increased until ∼2002 and then leveled off for the duration of the study. This increase probably was a response to the completion of major restoration activities in 1999 and corresponding improvements in water quality (Fig. 1). SiZer analysis identified a significant threshold response for abundance of metal-sensitive taxa in 2002, ∼10 y after restoration activities were initiated. Mean abundance of sensitive taxa at station AR3 approached values observed at the upstream reference station, a result suggesting that recovery had occurred for this particular metric. Additional remediation at this site probably will provide progressively smaller ecological benefits as the system continues to recover.

In addition to these long-term patterns, biological responses to improvements in water quality in the Arkansas River were influenced by both seasonal and annual variation in hydrologic characteristics. Although abundance of metal-sensitive species increased over time, recovery was briefly interrupted in the mid-1990s. This short-term temporal variation was very obvious in the multivariate data, which essentially tracked changes in the entire benthic community and showed several distinct thresholds (Fig. 3A–C). The abrupt shifts in community composition coincided with significant increases in annual stream discharge and heavy metals loading from 1995–1997. These changes were associated with several of the heaviest snow packs observed in our 17-y study ( http://www.wcc.nrcs.usda.gov/snow/) and highlight the importance of assessing restoration success within the context of episodic events and long-term variation in regional climate.

The location of an ecological threshold identified by SiZer is influenced by h, which essentially reflects the fit of the data to the LWPR model. In the example showing abundance of metal-sensitive taxa, the recovery threshold would occur ∼2 y later at h  =  2.5 (log₁₀h  =  0.4) and no thresholds would be observed at larger values of h (log₁₀h > 0.5). We suggest that a no-threshold model is unrealistic for these data because it does not reflect the obvious transition in 2002 and instead shows a monotonically increasing function where abundance eventually would exceed reference station values. Similarly, a larger h for the multivariate data (e.g., log₁₀h > 0.4) also would suggest a monotonically increasing function over time, which clearly does not reflect the complexity of this relationship. These different values of h effectively define different levels of temporal resolution of our data. With larger values we see long-term decadal patterns, but the annual and multiyear patterns are smoothed because the LWPR function is using data combined across many years (as indicated by the effective width of the LWPR). For smaller values of h (h  =  1.5), these multiyear patterns become more apparent, but annual variation is smoothed. Last, unrealistically small values of h (e.g., <1.0) would over-fit the data and produce a function that primarily tracks annual variability.

Water-quality criteria for lotic ecosystems are generally constructed as threshold values, and an exceedance of the threshold triggers an evaluation of sources and restoration alternatives (Davis and Simon 1995, Smith et al. 2005). However, this single-threshold approach to regulation might be inadequate for streams simultaneously exposed to natural (e.g., hydrologic variation) and anthropogenic (e.g., metals, sediment, and nutrients) stressors, where multiple thresholds can occur across time and space (Poole et al. 2004). Threshold models, such as piecewise linear analysis (Dodds et al. 2010), will always identify a single threshold, regardless of the nature of the stressor–response relationship or the number of thresholds that actually exist. In contrast, SiZer easily handles multiple thresholds and can provide important insights into the general form of stressor–response relationships or recovery trajectories.

Because of the inherent uncertainty in predicting threshold concentrations for contaminants, a more conservative approach is to identify a range of stressor values where threshold responses are likely to occur (Muradian 2001). We believe that SiZer is especially appropriate in situations where incorrectly identifying a threshold has serious consequences for protecting aquatic resources (e.g., overestimating a critical threshold concentration for a contaminant). SiZer displays thresholds as a function of various LWPR bandwidths and allows an investigator to consider a range of thresholds that reflect statistical uncertainty in the data at different levels of resolution. The SiZer map provides a transparent graphical representation of the threshold and how it responds to variability of the data. The SiZer approach also will be valuable in ranking sensitivity of ecological endpoints typically used for regulation and restoration goals. For example, compared to the patterns of community composition based on the multivariate data, abundance of sensitive species was considerably more variable and showed a smaller response to the episodic increases in Zn observed in the mid 1990s. These results suggest that community composition was more sensitive than abundances of sensitive species to short-term variation in stressor levels and might represent a more effective endpoint for assessing restoration success.