Defining ecosystem thresholds for human activities and environmental pressures in the California Current

Abstract

The oceans are changing more rapidly than ever before. Unprecedented climatic variability is interacting with unmistakable long‐term trends, all against a backdrop of intensifying human activities. What remains unclear, however, is how to evaluate whether conditions have changed sufficiently to provoke major responses of species, habitats, and communities. We developed a framework based on multimodel inference to define ecosystem‐based thresholds for human and environmental pressures in the California Current marine ecosystem. To demonstrate how to apply the framework, we explored two decades of data using gradient forest and generalized additive model analyses, screening for nonlinearities and potential threshold responses of ecosystem states (n = 9) across environmental (n = 6) and human (n = 10) pressures. These analyses identified the existence of threshold responses of five ecosystem states to four environmental and two human pressures. Both methods agreed on threshold relationships in two cases: (1) the winter copepod anomaly and habitat modification, and (2) sea lion pup production and the summer mode of the Pacific Decadal Oscillation (PDO). Considered collectively, however, these alternative analytical approaches imply that as many as five of the nine ecosystem states may exhibit threshold changes in response to negative PDO values in the summer (copepods, scavengers, groundfish, and marine mammals). This result is consistent with the idea that the influence of the PDO extends across multiple trophic levels, but extends current knowledge by defining the nonlinear nature of these responses. This research provides a new way to interpret changes in the intensities of human and environmental pressures as they relate to the ecological integrity of the California Current ecosystem. These insights can be used to make more informed assessments of when and under what conditions intervention, preparation, and mitigation may enhance progress toward ecosystem‐based management goals.

Introduction

Any performance management system—whether it is focused on health care, national security, business, or the environment—rests first and foremost upon a reliable set of measurements. These measurements, or indicators, can be used to evaluate a system's status (Otley ). The rise of ecosystem‐based management (EBM) of the ocean is no exception to this generality. As scientists and managers have tried to move marine EBM from theory to practice (Pikitch et al. , Arkema et al. , Levin et al. , Samhouri et al. , Link and Browman ), a frequent first step has been to identify and report on indicators of ocean conditions (Kershner et al. , Halpern et al. , Andrews et al. ). In marine EBM and elsewhere, however, indicators are necessary but not sufficient for successful performance management (Otley ). Useful indicators must describe a system in relation to the objectives set for it (Melnyk et al. ). Reference points provide just such a basis for comparison of an indicator's value to an internal or external standard. (See Appendix S1: Table S1 for definition of indicator, reference point, and other terms used in the text.) Thus, they are a critical component of advancing EBM from concept to commonplace.

An increasing variety of environmental performance management approaches are built around reference points, from global to local scales (Halpern et al. , Hamel et al. , UMCES ). For example, the most familiar reference point in fisheries science is the concept of maximum sustainable yield, which is calculated based on estimates of population size, carrying capacity, and growth rates (Hilborn and Walters ). One of the most informative techniques to establish reference points relies upon thresholds in ecosystem state. Following Groffman et al. (), we define an ecosystem threshold as a large, nonlinear change in an ecosystem state indicator (e.g., biodiversity of a key species or functional group) in response to an incremental change in an anthropogenic or environmental pressure(s), such as pollution or temperature (Lackey , Methratta and Link , Martin et al. , Samhouri et al. ).

Determining where ecosystem thresholds occur, and how large a shift may be induced by crossing them, is key for informing management that prepares for and is designed to avoid abrupt and undesired changes (Doak et al. ). For instance, the concept of rising variance, where the variability in an ecosystem state(s) can provide advance warning of an ecosystem shift, has led to new insights about both terrestrial and marine systems (Carpenter and Brock , Sydeman et al. , Litzow and Hunsicker ). Many analytical techniques for identifying ecosystem thresholds based on state–pressure relationships have also been proposed (Samhouri et al. , Large et al. , , b, Hunsicker et al. ), leaving an open question: Which is most appropriate, and under what conditions?

Here, we propose a framework centered on multimodel inference (MMI), rather than a single statistical tool, to define ecosystem thresholds for environmental and human pressures. To illustrate how our MMI framework can be used, we apply it to the U.S. California Current System (CCS), which supports more than $23 billion in revenue from fisheries, tourism, and recreation (data source: NOAA Coastal Services Center, 2013 GDP data for living resources and tourism and recreation sectors; http://www.oceaneconomics.org/). The vast majority of ocean condition reports for the CCS focus on ecosystem states, as well as environmental and human pressures, relative to internal standards—long‐term averages of conditions (Halpern et al. , Harvey et al. , Leising et al. , ONMS , Zador ). Here, we seek to enhance these assessments by directly connecting variability in environmental and human pressures to the potential for nonlinear ecosystem responses. Such nonlinear ecosystem responses are common in many ecosystems (Hunsicker et al. ) and can alter the human‐derived benefits of a system (Selkoe et al. ). As such, we believe our MMI framework that identifies ecosystem thresholds has broad application to other regions in the process of operationalizing EBM.

Methods

Analytical framework

We developed an analytical approach to define ecosystem thresholds for environmental and human pressures. The goal was to represent the value of a pressure relative to an inflection point in its relationship to one or more ecosystem states, and to quantify the magnitude of ecosystem change associated with crossing that value. This pressure–state approach is distinct from one that focuses on changes in state or pressure variables over time, or the time at which a threshold was crossed (Bestelmeyer et al. , Scheffer et al. ). The workflow (Fig. ) can be summarized in four parts, including (1) pre‐treatment, (2) screening, (3) functional form identification, and (4) threshold identification.

Analytical framework for defining ecosystem‐based thresholds for environmental and human pressures. S = ecosystem state indicator(s); E = environmental pressure indicator(s); H = human pressure indicator(s); DFA = dynamic factor analysis; mag = magnitude of ecosystem response across a threshold. Note that the tools listed here are intended as examples, rather than an exhaustive list. GAM, generalized additive model.

Step 1: Pre‐treatment: Which data?

The first step in this workflow represents scoping, winnowing, and data preparation. Scoping precisely defines the focus of the analysis in terms of ecosystem state and pressure indicators. It is important to clarify whether the focus of the analysis will be on univariate or multivariate ecosystem indicators, and on bivariate pressure–state relationships or associations between univariate ecosystem indicators and multiple pressures.

The goal of pre‐treatment is to narrow the universe of possibilities to a manageable subset of time series, a suite of spatial data from a limited time period, or a spatio‐temporal data set. This winnowing of data sources can be accomplished via expert opinion, reference to pre‐established indicator frameworks (e.g., European Marine Strategy Framework Directive, Tett et al. , California Current Integrated Ecosystem Assessment, Harvey et al. , the Puget Sound Partnership Science Update, PSP ), or multivariate analyses intended to reduce dimensionality for large data sets (e.g., dynamic factor analysis; Zuur et al. ). This step helps to reduce the occurrence of statistically spurious results. The process of winnowing also provides an opportunity to ensure that the state and pressure variable data are derived at spatial scales that allow for biologically plausible relationships. While it is difficult to provide prescriptive guidelines about requirements without a priori knowledge of effect sizes, variances, and other characteristics of the data, a minimum of ten matched state–pressure data points is recommended.

After winnowing, any data preparation such as the interpolation or extrapolation of missing data points should be completed. (We caution against using interpolated missing values when they occur at the beginning or end of time series). Decisions about whether and how to consider time lags and spatial‐scale mismatches in pressure–state relationships should be made at this stage. For example, are large‐scale indicators of the intensity of human activities appropriate for identifying threshold changes in an ecosystem state measured at a smaller spatial scale? Similarly, are thresholds in a pressure–state relationship expected to occur simultaneously, or would thresholds be more likely to be observed when ecosystem state data are lagged?

Step 2: Screening: Are there nonlinearities?

The second step explores the potential for nonlinear relationships between ecosystem states and pressures and seeks to identify the potential existence of thresholds. This portion of the analysis relies on MMI, rather than a single statistical tool. Broadly defined, MMI is the application of several quantitative representations of a system to learn how the system works (Townsend et al. ). As this step is intended to eliminate pressure–state relationships from the analysis that are unlikely to show evidence for thresholds, we recommend comparing model results qualitatively (rather than developing a quantitative model consensus). Strongest inferences can be made where models agree on the existence and location of thresholds, and quantitative descriptors of the identity and magnitude of the threshold can be developed in the following steps.

Common approaches for identifying the presence of thresholds include gradient forest, generalized additive and generalized additive mixed models (GAMs and GAMMs), and specified functional form analyses (Samhouri et al. , Large et al. , , b, Baker and Hollowed , Hunsicker et al. ). Importantly, the choice of methods may be influenced by the completeness of the indicator data, as some methods use all ecosystem and pressure indicators simultaneously and require complete data sets (e.g., gradient forest), while other methods evaluate individual pressure–state relationships and can reasonably handle missing values in a time series (e.g., GAM).

Step 3: Functional form identification: What type of nonlinearity exists?

For relationships identified as nonlinear in the preceding screening step, the next stage of analysis derives relevant statistics from the models to describe their signs and functional forms. The outcome of this analysis should be a quantitative description of whether a pressure is positively or negatively correlated with an ecosystem state, and at a minimum a qualitative description of the shape of the relationship (e.g., hockey stick, sigmoidal, parabolic). While many statistical models will allow such a quantitative description (e.g., GAMs, specified functional forms), others will not (gradient forest analysis).

Step 4: Threshold identification: How strong are the nonlinearities?

For relationships that emerge as nonlinear, a final set of analyses is used to determine the location (inflection point) and strength of the threshold. This step provides (1) a quantitative estimate of the threshold level(s) of a pressure corresponding to an abrupt change in the direction of its relationship with an ecosystem state, typically defined as the point of inflection where the second derivative changes sign (Samhouri et al. , Large et al. ), and (2) the magnitude of change in an ecosystem state associated with breaching the threshold level of a pressure. Several statistical tools can be used to locate the threshold, estimate the uncertainty around its location, and describe the magnitude or effect size corresponding to the threshold (Andersen et al. , Foley et al. ).

California current application

We applied the MMI framework described above to the CCS. A comprehensive assessment of the conditions within this coupled social–ecological system is reported in NOAA's California Current Integrated Ecosystem Assessment (CCIEA; Harvey et al. ). The following section represents the pre‐treatment steps for application of our MMI framework for identifying thresholds.

The CCS is an eastern boundary current ecosystem, with seasonal periods of upwelling of cold nutrient‐rich waters along the coast that drive primary and secondary productivity and affect the dynamics of the diverse resident and migratory species throughout the food web (Bograd et al. , Hazen et al. ). In addition to fluctuations in the oceanographic and biophysical environment, the CCS is affected by a variety of human uses that include fisheries and other marine activities, as well as land‐based activities that result in localized (i.e., point source) and broadscale (i.e., nonpoint source) transfer of materials to the coastal zone (Halpern et al. , Andrews et al. ).

Changes in environmental pressures in the CCS can be abrupt. Some vary at relatively short timescales, including short‐term variability in upwelling strength as tracked by the North Pacific High (Schroeder et al. ) and the Northern Oscillation Index (NOI; Schwing et al. ). Other processes act at interdecadal scales; for example, the strength of transport by the North Pacific Gyre is indexed by the North Pacific Gyre Oscillation (NPGO; Di Lorenzo et al. ), and decadal changes in sea surface temperature regimes in the northeast Pacific are tracked by the Pacific Decadal Oscillation (PDO; Mantua and Hare ). When warmer‐than‐average temperatures and weak upwelling dominate the CCS (e.g., positive PDO and negative NOI), large ecosystem state changes have been observed (King et al. ), including shifts in planktonic communities (Peterson et al. ) and lower trophic‐level fishes (Chavez et al. ), as well as higher trophic‐level fishes (Lindley et al. ), seabirds, and mammals (Leising et al. ).

Human activities within the CCS are also quite dynamic (Andrews et al. ). Economic shocks (e.g., McKenna et al. ), emerging technologies (Kim et al. , Plummer and Feist ), and regulatory shifts (e.g., in fisheries; Hilborn et al. , Lubchenco et al. ) have caused rapid changes in ocean uses. As with variability in the environment, these changes have both direct and indirect effects on various components of the ecosystem (Halpern et al. ). Here, we drew from the 19 human activities presented in the CCIEA (Harvey et al. ) and focused on the 10 that provided available data across most of the period of interest (Table ).

Indicators of environmental and human pressures in the California Current System

Component	Pressure	Indicator	Time series
Environmental	Basin‐scale sea surface temperature	PDO summer and winter indices	1900–2014
Basin‐scale atmospheric forcing	NOI summer and winter indices	1948–2014
Changes in source waters	NPGO summer and winter indices	1950–2014
Human	Atmospheric pollution	Deposition of sulfate	1994–2014
Commercial shipping activity	Volume of water disturbed	2001–2012
Dredging	Dredge volumes	1997–2014
Groundfish fisheries removals	Commercial groundfish landings	1981–2014
Habitat modification	Distance trawled	1999–2012
Inorganic pollution	Toxicity‐weighted chemical releases	1988–2012
Invasive species	Tons of cargo	1993–2013
Nutrient input	Nitrogen and phosphorus input	1945–2010
Organic pollution	Toxicity‐weighted concentrations of pesticides	1993–2010
Total fisheries removals	Total commercial and recreational landings	1981–2014

Note

Further details of data sets can be found in Appendix S1: Table S2.

Data sources

Our analysis centered on the identification of levels of pressures likely to induce the crossing of a threshold for one or more indicators of ecological integrity and was based on time series data. While spatial variability is certainly a key feature of the California Current, an analysis of whether and how relationships between ecosystem states and pressures differ from place to place was beyond the scope of this study.

The CCIEA defines ecological integrity as the species composition, biodiversity, and functional organization of an ecosystem and includes nine coast‐wide indicators: a mean trophic index, species density, species richness, Simpson diversity, scavenger biomass, the Northern copepod anomaly (winter and summer), and California sea lion (Zalophus californianus) pup abundance and pup growth rate (Table ; Williams et al. ). These ecosystem state time series were derived from three fishery‐independent data sets that span several major taxonomic groups (invertebrates, fishes, and mammals), and summarized at the largest spatial domain possible given the data set (see Williams et al. for details).

Ecosystem states assessed in the case study of the California Current System

Ecosystem attribute	Indicator	Taxonomic group	Community	Definition and source of data	Time series	Sampling frequency
Species composition	Northern copepod anomaly, winter and summer	Invertebrates	Pelagic	Monthly anomalies in the relative biomass of copepods with cold‐water affinities off Newport, OR (Peterson et al. , NOAA)	1996–2014	Biweekly; summarized as winter and summer anomalies
Biodiversity	Species richness, Species density, Simpson's index	Groundfish	Benthic	Index of groundfish community composition (Bradburn et al. , NWFSC)	2003–2014	Summers, annual
Functional organization	Mean trophic index	Groundfish	Benthic	Trophic structure of groundfish community (Bradburn et al. , NWFSC)	2003–2014	Summers, annual
Functional organization	Scavenger biomass	Groundfish and invertebrates	Benthic	Relative biomass of scavengers, as defined by esp. Brand et al. (), from fishery‐independent surveys (Bradburn et al. , NWFSC)	2003–2014	Summers, annual
Functional organization	California sea lion pup production and growth	Marine mammals	Top predators	Average no. of pups on San Miguel Island in late July and predicted daily growth rate of pups between June and October (Melin et al. , Wells et al. , NMML)	1997–2014	July and June–October, respectively; annual

Notes

These “Ecological integrity” indicators are used by NOAA's Integrated Ecosystem Assessment for the California Current with detailed descriptions in Williams et al. (). Note that the California sea lion pup growth time series was missing data for a single year (2011). Preliminary investigation suggested that a simple mean interpolation was the most parsimonious way to fill this gap. NOAA = National Oceanic and Atmospheric Administration; NMML = National Marine Mammal Laboratory; NWFSC = Northwest Fisheries Science Center.

There are a wide variety of pressures with the potential to effect change in these ecosystem states. Although there are regional differences within the CCS in physical forcing, climatic variability, and ecosystem responses (Mendelssohn et al. , García‐Reyes and Largier ), we focused on six basin‐scale environmental pressures in the North Pacific (Table ; Appendix S1: Table S2). These environmental pressures are the primary basin‐wide indicators of oceanography and climate in the CCIEA (Hazen et al. ) and include the winter and summer anomalies in the NOI (NOI_w and NOI_s, respectively), NPGO (NPGO_w and NPGO_s, respectively), and PDO (PDO_w and PDO_s, respectively). We include both winter and summer indicators because two modes of upwelling in the CCS drive different components of biological productivity. While the strongest upwelling occurs in the summertime, the winter mode (upwelling at the beginning of the season) can be equally important for growth and reproduction in some species (Black et al. ).

As with the environmental pressures, many human activities are likely to be associated with variability in ecological integrity in the CCS. Here, we focus on ten human activities that have received global (Halpern et al. , ) and regional (Halpern et al. , Andrews et al. ) attention in relation to changes in marine ecological integrity (Table ; Appendix S1: Table S2). Four activities relate to pollution (atmospheric, inorganic, nutrients, and organic), three pertain to habitat disturbance (commercial shipping, habitat modification, and dredging), two concern extraction (total fisheries and groundfish landings), and one is associated with invasive species. The intensities of all of these activities were summarized at the scale of the full CCS (see Andrews et al. for details).

Implementation of the analytical framework

To embrace the MMI philosophy of the framework described above, we used gradient forest, GAM, and GAMM analyses to screen for nonlinearities and potential threshold responses of ecosystem states (Liaw and Wiener , Ellis et al. , Baker and Hollowed , Large et al. , Hunsicker et al. ). We tested for nonlinearities in all possible pressure–state relationships, but subsequently excluded those without plausible mechanistic relationships. While this approach increased the possibility that we would detect significant thresholds and nonlinearities, our interest was in the precautionary identification of potential thresholds rather than statistical significance (White et al. ), especially given the limitations of inferences based on P‐values (Barber and Ogle ).

We implemented the gradient forest, GAM, and GAMM analyses on a subset of the full time series. When considered collectively, the time series of ecosystem states spanned a 19‐yr period from 1996 to 2014 (Table ), but the span of overlap among all ecosystem and pressure time series was 10 yr (2003–2012). For the gradient forest analyses, missing data for any time series is problematic; we truncated each time series to this 10‐yr period for the purposes of MMI (hereafter, truncated analyses). We removed two human pressures (organic pollution and nutrient input) from the truncated analyses due to missing data in 2011–2012.

While the truncated analyses allowed for direct comparison of results among models, this approach had the obvious drawback of not using all the available data. Therefore, we also conducted GAM/GAMM analyses on the longest time series available for each pressure–state pair (hereafter, full GAM analyses).

Detailed description of the gradient forest analyses

The gradient forest analysis quantified the ability of each environmental or human pressure to predict variation in the time series of multiple ecosystem states (Breiman , Large et al. ). Gradient forest analysis is an ensemble of random forest models, each of which splits ecosystem states into two groups at specific values of an environmental or human pressure. Partitions are further made until one group becomes a terminal node. The R²‐importance of each value of the pressure—in other words, the possibility that it represented a threshold—was calculated based on maximizing the homogeneity of variance of the ecosystem state values within each subsequent partition. Ecosystem states with zero variance explained by random forests were not included in the final gradient forest models.

We estimated the aggregate response of all ecosystem indicators, called the cumulative ecosystem response, to each pressure. The cumulative ecosystem response was calculated from the cumulative importance distributions of split improvement for each ecosystem state, scaled for each ecosystem indicator by R²‐importance and standardized by the density of observations. Though the gradient forest analysis identified potential thresholds wherever splits in the pressure improved the homogeneity of variance of the cumulative ecosystem response values, we focused only on improvements in cumulative ecosystem response R²‐importance that exceeded 0.01.

For pressures that increased cumulative ecosystem response R²‐importance by ≥0.01, we identified individual ecosystem state thresholds based on a range of pressure values. This range was defined in relation to significant increases in the cumulative R²‐importance values of each ecosystem state given a split at specific values of the pressure (Large et al. ). We used the R packages “randomForest” (Liaw and Wiener , R Core Team ) and “gradientForest” (Ellis et al. ) for all calculations.

Detailed description of the GAM/GAMM analyses

We used a model selection approach to determine whether a nonlinear GAM or GAMM provided a more parsimonious explanation of pressure–state relationships than a linear model (Sonderegger et al. , Bestelmeyer et al. , Samhouri et al. ). However, because our analyses focused on time series data, we first used a log‐likelihood ratio test adapted from the methods described in Gilman et al. () to determine whether a generalized additive mixed model with autocorrelated error structure (GAMM) was more appropriate than a GAM with normal error structure. The log‐likelihood ratio test was based on a comparison of (1) the fit of residuals from the GAMM in a linear model with autocorrelated structure in the residual covariance matrix to (2) the fit of residuals in a simple linear model. Where temporal autocorrelation was determined to be significant, we used Akaike's Information Criterion (corrected for small sample sizes; AIC_c) to select between GAMM (“gamm” function in R) and linear models with autocorrelated structure in the residual covariance matrices (Eqs. ). Where temporal autocorrelation was deemed nonsignificant, we used AIC_c to select between simple GAM (gam() function in the mgcv() package in R) and linear models (Eqs. , and ).

Formally, let E_y be the value of the ecological indicator in year y, α a fixed intercept effect, P_y the value of the pressure in year y, s() the smoothing function, ε a normally distributed random error term, and R_ac and R the structures of the residual covariance matrices (represented here for a three‐year period), with and without autocorrelation, respectively, yielding:[Image Omitted. See PDF][Image Omitted. See PDF][Image Omitted. See PDF][Image Omitted. See PDF]

The formulation of Eq.  implies that the most recent values of a pressure have the greatest influence on the response of the ecosystem indicator, and the previous influence of a pressure diminishes quickly according to the autocorrelation coefficient ρ. In contrast, the formulation of Eq.  implies that values of a pressure in a previous year have zero influence on the response of an ecosystem indicator.

We determined whether there was more evidence for linearity or nonlinearity based on the following three criteria: (1) The estimated degrees of freedom, which increases with a curve's nonlinearity, was ≥2.0 (Zuur et al. , Hunsicker et al. ); (2) the generalized cross validation score (GCV; Wood ) was minimized in the GAM compared to the LM (GCV is not a calculated criterion for GAMMs); and (3) the difference in AIC_c values was ≥2.0 units in favor of the nonlinear model (Burnham and Anderson ). For both GAM and GAMM, we used thin plate regression spline smoothing terms, the “mgcv” (Wood ) package, and set the size of the basis dimension to 4 to reduce the possibility of over‐fitting (Large et al. ).

For pressure–state relationships identified as nonlinear, we defined the location of the threshold as the inflection point, that is, the value of the pressure where the second derivative changed sign (Fewster et al. , Bestelmeyer et al. , Samhouri et al. , Large et al. ). For these analyses, we calculated the 95% CI of the smoothing function itself, along with its second derivative, via bootstrapping of the residuals in order to allow for autocorrelation. This procedure generated a range of pressure values where a threshold might occur. Because location of the inflection point on the second derivative function is based on a statistical fit and cannot be determined exactly, we defined that location as the place where the second derivative is most different from zero. Thus, we are most confident that the second derivative changed from zero to nonzero at this point.

We determined the functional form and magnitude of change in ecosystem state variables associated with crossing a threshold based on the full GAM/GAMM analyses only. The functional form of the nonlinear GAM/GAMMs was defined based on the sign of the smoother coefficient on each pressure, the number of inflection points, and visual inspection of the shape of each pressure–state relationship. The magnitude of change associated with crossing a threshold was estimated based on the proportional difference in ecosystem indicator values on either side of the threshold, bounded as the maximum, minimum, and best estimate of the threshold value (as defined in the preceding paragraph).

Model output

In the interest of using MMI to detect thresholds, we report the results of each of the three analyses (gradient forest analysis, full GAM analyses, and truncated GAM analyses). Note that for both the full and truncated GAM analyses, model selection criteria never identified GAMs with autocorrelated error structure (GAMM) as the most parsimonious model; thus, the following results are based on output from GAMs.

For relationships with identified thresholds and plausible mechanistic explanations, we mapped the range of threshold values for each ecosystem indicator to the corresponding pressure time series to identify the relative frequency with which the pressure exceeded the ecosystem‐based threshold.

Results

Our analysis reveals effects of both environmental and human pressures as major predictors of individual and aggregated ecosystem states in the CCS. These threshold relationships were evident for individual ecosystem indicators as well as the cumulative ecosystem response, though the types and magnitudes of nonlinearities varied among ecosystem states. Encouragingly, while the importance of many of these pressures for the CCS was previously known, we identified several novel nonlinear pressure–state relationships with potential utility for EBM. Below we present the results such that they parallel the MMI framework introduced in Fig.  (except for pre‐treatment steps described in ). We also demonstrate how the ecosystem thresholds can be mapped to time series of pressures, providing new insights into how variability in the environment and human activities may influence ecological integrity.

Screening for nonlinearities

We identified two pressure–state relationships as potentially nonlinear via multiple methods. Both the truncated and full GAM analyses suggested that the relationships between (1) the winter copepod anomaly and habitat modification, and (2) sea lion pup production and PDO_s, were characterized by thresholds. In addition, we flagged nine other potential nonlinear pressure–state relationships based on a single method (Figs. , Table ; Appendix S1: Table S3).

Gradient forest analysis of nine ecosystem states, six environmental pressures, and eight human pressures using data from 2003 to 2012. Top: importance of environmental and human pressures weighted across all ecosystem states. Pressures with R2‐importance ≥0.01 were considered capable of predicting variation in ecosystem states. Cumulative importance (in R2‐importance units) of the aggregate response of all ecosystem states (middle) and four individual ecosystem states that were predicted by the best model for this set of pressures (bottom) across the gradient of each pressure. Each plot is scaled to the maximum cumulative response to allow for direct comparison of ecosystem responses to each pressure. PDO, Pacific Decadal Oscillation; NOI, Northern Oscillation Index; NGPO, North Pacific Gyre Oscillation; N copepod anom_s, Northern copepod anomaly summer; N copepod anom_w, Northern copepod anomaly winter.

Truncated GAM analyses (data from 2003 to 2012) of ecosystem responses to environmental or human pressures, where the dashed line is the GAM smoother, gray polygon is 95% CI, points are raw data, thick solid line indicates the threshold range where the 95% CI of the second derivative does not include 0, and red dotted arrow indicates the best estimate of the location of the threshold (i.e., where the second derivative is most difference from zero within the threshold range). See Appendix S1: Fig. S4 for additional details.

Full GAM analyses (data from 1996 to 2014) of ecosystem responses to environmental or human pressures. See Fig.  caption and Appendix S1: Fig. S5 for additional details.

Nonlinear ecosystem state–pressure relationships as identified by gradient forest (GF), truncated generalized additive model (GAM), or full GAM analyses including functional form, location of threshold range, and magnitude of change in ecosystem state across the threshold range

Ecosystem state	Driver/Pressure	Analysis	Functional form(s)	Location of threshold(s)	Best estimate of threshold location(s)	Magnitude of response(s) (%)
Copepod anomaly winter	PDO winter	GF	–	−0.5 to −0.2	–	–
Copepod anomaly winter	Habitat modification	Truncated GAM	Parabolic	143–234	208	70
Full GAM	Sinusoidal	138–252	227	30
Copepod anomaly summer	NPGO winter	GAM	Hockey stick	0.2–0.8	0.2	180
Copepod anomaly summer	PDO summer	GF	–	−1.2 to 0.5	–	–
Copepod anomaly summer	PDO winter	GF	–	0.7–0.8	–	–
Scavenger ratio	Commercial shipping activity	GF	–	14.7–15.2	–	–
Scavenger ratio	PDO summer	GF	–	−0.6 to 0.1	–	–
Groundfish mean trophic index	PDO summer	GF	–	−0.3 to 0	–	–
CA sea lion pup production	NOI summer	GAM	Hockey stick	−0.4 to 1.2	0.2	10
CA sea lion pup production	PDO summer	Truncated GAM	Sigmoidal	−1.5 to −0.2	−0.8	10
Full GAM	Hockey stick	NTI	NTI	NTI
CA sea lion pup production	PDO winter	Truncated GAM	Sigmoidal	0.7–1.5	0.9	30
Truncated GAM	−1.4 to 0.2	−0.8	0

Note

NTI, no threshold identified by CI of the second derivative; –, information not determined from the model.

PDO_s and PDO_w both showed evidence for nonlinear relationships with three ecosystem states and the cumulative ecosystem response, underlining the potential importance of environmental conditions for the CCS. In contrast, NPGO_w, NOI_s, and habitat modification each only showed evidence for nonlinear relationships with a single ecosystem state. While both gradient forest and GAM analyses pointed to individual ecosystem thresholds in relation to PDO_s and PDO_w, only gradient forest analysis highlighted a potential ecosystem threshold in relation to commercial shipping activity (Fig. , Table ; Appendix S1: Tables S3 and S4).

There was no evidence for nonlinear responses of groundfish species density, species richness, or Simpson diversity to the environmental and human pressures we tested. Of the five ecosystem states with nonlinear responses, four showed evidence for thresholds in relation to more than one pressure (Table ; Appendix S1: Table S3).

Functional form identification

The truncated and full GAM analyses allowed identification of functional forms of relationships between three ecosystem states and five pressures. California sea lion pup production showed evidence for hockey stick and inverse parabolic relationships with three environmental pressures. First, it declined precipitously with initial increases in PDO_s but thereafter was relatively invariant (Figs.  and ). This ecosystem state also showed an inverse parabolic relationship with PDO_w, increasing initially to a peak at a value of around −0.5, and then declining with large positive values of this pressure (Fig. ). Finally, California sea lion pup production exhibited little change associated with increases in NOI_s until NOI_s values of 0.25, whereby increasing NOI_s values were associated with anomalously low pup production rates (Fig. ).

We also determined the functional form of nonlinear relationships between the summer and winter copepod anomalies and both environmental and human pressures. The relationship between the summer copepod anomaly and the NPGO_w was best described as a hockey stick, such that this ecosystem state was associated positively and linearly with NPGO_w when NPGO_w was negative, but was relatively invariant for positive values of NPGO_w (Fig. ). The winter copepod anomaly showed a parabolic relationship with habitat modification, with maximum values at low and high values of this pressure, according to the truncated GAM analysis (Fig. ). However, the full GAM analysis indicated precipitous declines in the winter copepod anomaly associated with initial increases in habitat modification, followed by small increases at intermediate values of habitat modification, and finally a second set of declines associated with high levels of habitat modification (similar to a sinusoidal relationship; Fig. ). The contrast in the functional form of this pressure–state relationship demonstrates the potential for strong leverage of a small number of extreme data points when time series are relatively short.

Quantifying the strength of nonlinearities

We quantified the location of thresholds that emerged from both the gradient forest and GAM analyses, and the magnitude of change in ecosystem states associated with breaching threshold levels of pressures based on the GAM analyses. The values of environmental pressures associated with threshold changes in ecosystem states corresponded to changes from positive to negative anomalies in some cases (winter copepod anomaly and California sea lion pup production vs. PDO_w; summer copepod anomaly, scavenger ratio, and groundfish mean trophic index vs. PDO_s; California sea lion pup production vs. NOI_s), but not in all comparisons (summer copepod anomaly vs. NPGO_w and PDO_w; California sea lion pup production vs. PDO_s; Table , Fig. ; Appendix S1: Figs. S4 and S5). One nonlinear relationship did not meet our definition of having a distinct threshold (CA sea lion pup production vs. PDO_s full GAM analysis).

Location of thresholds for multiple ecosystem indicators related to (a) the summer mode and (b) the winter mode of the Pacific Decadal Oscillation (PDO). Open circles: gradient forest; closed circles: GAM. Note that for CA sea lion pup production in relation to the PDO winter, there are two key thresholds at both cold and warm anomalies (also see Table ).

For the thresholds related to human pressures (Table , Fig. ; Appendix S1: Figs. S4 and S5), the truncated and full GAM analyses agreed on the location of the threshold for the winter copepod anomaly in relation to habitat modification. The threshold for commercial shipping activity determined by the gradient forest analysis occurred at intermediate values of this activity.

We observed a range of magnitudes in the response of ecosystem states across the associated thresholds. At the low end, California sea lion pup production changed by ~10% across the threshold found with NOI_s (Fig. , Table ). At the high end, the winter northern copepod anomaly changed by 180% across the threshold found with NPGO_w (Fig. , Table ).

Mapping ecosystem threshold values onto pressure time series

A useful way to visualize how the temporal changes in the California Current pressures relate to ecosystem states is to plot the values of pressures relative to ecosystem thresholds (Fig. ). This approach complements one that relies strictly on anomalous variation to determine whether conditions in any particular year are “good” or “bad.” We focus here on threshold values of pressures based on the full GAM analyses, which harnesses all of the data available to us. For example, in 10 of 19 yr between 1996 and 2014 the NPGO_w suggested conditions distant from the threshold corresponding to an abrupt decline in the summer copepod anomaly (blue shading), one year in which the NPGO_w was in the range of the threshold for the summer copepod anomaly (yellow shading), and 6 yr in which the NPGO_w had exceeded this threshold (red shading; Fig. a). In nine of 19 yr between 1996 and 2014 the NOI_s index suggested conditions distant from the threshold corresponding to an abrupt decline in CA sea lion pup production (blue shading), one year in which the NOI_s was in the range of this ecosystem threshold, and one year in which the NOI_s exceeded the threshold corresponding to an abrupt decline in CA sea lion pup production (1998; Fig. b).

Time series of (a) the winter Northern Pacific Gyre Oscillation (NPGO) and (b) the summer Northern Oscillation Index (NOI), relative to confidence intervals for thresholds in ecosystem states (horizontal gray lines). The NPGO winter is shown in relation to the summer copepod anomaly, while the NOI summer time series is shown in relation to thresholds for California sea lion pup production (note that pup production is higher when values of NOI summer are lower). Blue shading indicates that the value of the environmental pressure was associated with more positive values of the ecosystem state, while red shading indicates the opposite. Yellow shading indicates that the value of the environmental pressure was within the confidence interval of the threshold for the relationship between it and the ecosystem state.

Discussion

The quantitative identification of abrupt changes in ecosystems is an essential step toward forecasting and preparing for the accelerating changes of the Anthropocene (Rockström et al. , Steffen et al. , Scheffer et al. ). In marine environments, the consequences of crossing thresholds—for population extinctions, anoxia, and acidification, for instance—are often difficult to reverse (deYoung et al. , Selkoe et al. ). Anticipating, mitigating for, or avoiding (when possible) these thresholds may reduce the risk of unwanted collapses and buffer coupled social–ecological systems from dramatic change (e.g., Groffman et al. ). Even in cases where environmental changes beyond the control of resource managers cause threshold ecosystem responses, such recognition provides crucial context for decisions about human uses of the marine environment that are subject to manipulation. The quantitative framework outlined here can help to define ecosystem‐based thresholds for human and environmental pressures, an issue with potentially broad application to other regions in the process of operationalizing EBM. Indeed, our application of the framework to the CCS revealed novel nonlinear ecosystem responses, consistent with well‐understood oceanographic forcing mechanisms.

Multimodel inference for ecosystem thresholds

Using indicators developed by NOAA's CCIEA, we screened all bivariate state–pressure relationships for nonlinearities using multiple analytical approaches. This strategy has clear strengths, as well as, arguably, some weaknesses. To its advantage, our framework introduces a common, stepwise, accessible, and data‐driven approach to characterizing thresholds in state–pressure relationships (cf. Bestelmeyer et al. ). As such, it can provide multiple lines of evidence for associations between certain pressures (e.g., PDO) and threshold responses in an array of ecosystem states. Our comprehensive exploration of thresholds also supports the concept that some ecosystem states are more prone to threshold responses than others (e.g., the northern copepod biomass anomaly; Fig. ). Furthermore, employing multiple statistical approaches helped to identify ecosystem thresholds that individual methodologies might have missed. As an example, only the full GAM detected a nonlinear, threshold response of California sea lion pup production to NOI_s (Fig. ). Though beyond the scope of the current study, we further suggest that this MMI approach could be used (1) to harness the power of spatio‐temporal data sets to test for threshold changes in ecosystem states across spatial gradients in pressures, and (2) to determine whether there are nonlinear relationships between ecosystem states, such as in the case of trophic relationships between the abundance of prey and their predators.

While our multimodel approach helped identify thresholds, using multiple lines of evidence sometimes led to inconsistent conclusions. For example, there was only one pressure–ecosystem state relationship for which thresholds were identified by multiple analyses (northern copepod anomaly in the winter with habitat modification—full and truncated GAM analyses; Table ). As mentioned above, a small number of extreme data points caused the functional form of this pressure–state relationship to differ in the analysis of the shorter and longer data sets. Furthermore, the mechanistic underpinnings for this relationship are unclear, though it is possible that indirect effects cause habitat modification (trawling) to lag changes in the northern copepod anomaly (ecosystem productivity). In addition, we tested >100 state–pressure relationships, increasing the possibility of making conclusions about thresholds that were statistically spurious. This choice was sensible because our interest was in the identification of potential thresholds rather than statistical significance (White et al. ), but may not be appropriate for all applications. We caution against using this exploratory approach to identify novel causal relationships. Instead, it is best used for generating new hypotheses or confirming existing ones regarding mechanistic links between ecosystem pressures and states. Lastly, we did not specifically examine threshold responses of multivariate ecosystem indicators or the potential for thresholds to emerge from multivariate predictors (Large et al. ), though we did identify threshold changes in the cumulative ecosystem response using gradient forest analysis (Fig. ). These avenues are ripe for future research, with multivariate thresholds providing context for evaluating tradeoffs. We caution that it may be challenging to use multivariate thresholds to inform management decisions (Fay et al. ).

Given these strengths and weaknesses, we believe the MMI approach may be particularly useful as a precautionary framework for identifying ecosystem components and relationships worthy of further detailed analyses. We contend that the risks of identifying spurious thresholds and over‐interpreting a limited data set are counterbalanced by the need to anticipate, prepare for, and, if necessary, avoid crossing ecosystem thresholds (sensu Jacquet et al. ).

Ecosystem thresholds in the California Current

In the CCS, we identified the existence of threshold responses of five ecosystem states to four environmental and two human pressures. These results broadly corroborate existing evidence for the importance of physical forcing in the CCS (Di Lorenzo et al. , Bograd et al. , Black et al. ). Previous studies have established the strong influence of human activities in this system (Halpern et al. , Hilborn et al. , Andrews et al. ). Our results add to this literature by demonstrating that while ecosystem states have changed in response to human activities over the relatively short length of our time series, there is little evidence that these changes have been nonlinear. However, this result does not imply that linear relationships between ecosystem states and human activities are inconsequential. Indeed, some may require increased attention by resource managers when environmental conditions are poor. For example, given that sea lion pup production declines precipitously when NPGO winter is negative, it may be important to provide increased scrutiny on human activities like fishing that affect sea lion prey during those periods (Chasco et al. ).

Nonetheless, the analyses presented here stand in contrast to results from other systems. For instance, in the Northwest Atlantic fisheries‐related pressures were associated with substantial threshold responses in the ecosystem (Large et al. , ). It is important to note that the Northwest Atlantic analysis focused specifically on ecosystem responses to fishing and that the lengths of the CCS time series we used were comparatively short. This latter point is a limiting factor in identifying ecosystem thresholds for many statistical methods and may explain why we did not observe nonlinear associations between ecosystem states and fisheries landings in the CCS.

Notwithstanding such limitations, we found two significant threshold relationships in the full GAM analyses with plausible mechanistic underpinnings. First, the summer northern copepod anomaly showed a positive, linear increase before asymptoting at values >0.2 with the winter mode of the NPGO (Fig. ). Increasingly positive values of the NPGO_w are associated with greater upwelling strength, nutrient transport into the photic zone, and primary productivity (Di Lorenzo et al. ), which could fuel the increased secondary productivity we observed. Second, California sea lion pup production showed a precipitous decline when values of the summer mode of the NOI exceeded 0.2 (Fig. ). This threshold response is counter to the assumption that highly positive NOI_s values (La Niña conditions) should be good for pup production. However, the NOI_s also captures the quick transition between El Niño and La Niña conditions experienced in the North Pacific, especially between 1998–99 and 2009–10. Thus, the 1998 and 2010 NOI_s values were the most positive values in the time series and represent the subsequent La Niña conditions of the following years. It was curious that the winter mode of the NOI did not identify this nonlinear response in the hypothesized direction. These relationships highlight the importance of considering time lags between the pressure and the ecological response, which should be the subject of further research.

Perhaps the most interesting outcome of applying the MMI framework to the California Current data is the strong evidence for threshold ecosystem responses to both modes of the PDO. Previous studies have identified climate‐driven “regime shifts” as recently as 1998, as the PDO reversed sign and many ecosystem components including salmon and anchovy responded (Peterson and Schwing ). While the importance of PDO in the California Current has a long history (Hare et al. , Mantua and Hare ), evidence for associated nonlinear ecosystem responses is more limited. For example, Bestelmeyer et al. () provided evidence for a linear relationship between euphausiid abundance and the PDO in the Southern California Current, and this conclusion has been complemented more recently by analyses showing that copepod species and coho salmon in the Northern California Current linearly track shifts in the PDO (Litzow and Hunsicker ). In contrast, our analyses imply that as many as five of the nine ecosystem states we evaluated exhibited threshold increases in response to negative PDO_s values (copepods, scavengers, groundfish, and marine mammals; Fig. a, Table ), a result consistent with the idea that the influence of the PDO extends across multiple trophic levels (Hare et al. , Mantua and Hare ). The location of thresholds associated with PDO_w for copepods and marine mammals ranged more widely (Fig. b, Table ), suggesting taxon‐specific responses. These interpretations are supported by results from the gradient forest analyses, which show that the aggregate ecosystem state in the CCS demonstrates a threshold relationship with both modes of the PDO (Fig. ).

Recognition of taxon‐specific vs. ecosystem‐wide threshold responses is a novel insight that can be tailored to sector‐specific or full‐ecosystem management needs. Instead of categorizing all anomalously positive (or negative) pressure values as “good” or “bad” across all ecosystem components, taxon‐specific thresholds allow a more refined evaluation of variability in environmental pressures (e.g., Fig. ). On the other hand, cross‐taxa agreement on the location of threshold values of a pressure provides insights into ecosystem‐wide responses to changing ocean conditions. In our application, gradient forest analysis identified a threshold cumulative ecosystem response to PDO_s (Fig. ) that spanned the same range of values as thresholds for individual ecosystem states identified via GAMs (Figs. ) and gradient forest analysis (Figs.  and ). These results clearly point to important nonlinearities in ecosystem dynamics across a range of trophic levels.

Conclusion

This paper outlined a quantitative framework based on MMI that allows for precautionary screening of threshold relationships between ecosystem states and environmental or human pressures. Establishment of quantitative ecosystem‐based thresholds allows for direct ecological interpretation of variability in environmental pressures and the intensity of human activities. This approach provides a foundation for evaluating the status of ecosystems and informing precautionary management of the multitude of human activities occurring within them. The most successful applications of this framework will provide advice to managers charged with developing ecosystem monitoring priorities, prepare ocean users for major shifts in ecosystem conditions, and indicate when human activity levels risk imposing critical transitions in individual or collective ecosystem states.

Applied to the CCS, we demonstrated how the MMI framework highlights potentially nonlinear ecosystem state–pressure relationships. The CCS is a system increasingly crowded with human uses (Halpern et al. ), yet our results suggest that while such activities may be associated with linear changes in the ecosystem over the time period of our study, broadscale oceanographic indices such as the PDO, the NPGO, and the NOI are primarily responsible for threshold changes. These findings provide support for emphasizing the importance of these particular environmental pressures to resource managers of the California Current. In particular, the PDO appears to be more associated with nonlinear ecosystem responses than the other pressures we examined. While resource managers cannot directly and easily control large‐scale environmental pressures, tracking their status and strength provides critical context for decisions about human activities that can be managed. Along with parallel analyses in other regions of the United States (Large et al. , ; Tam et al., unpublished manuscript), our results underscore the importance of long‐term monitoring data to capture ecosystem‐wide changes driven by environmental and anthropogenic pressures.

Acknowledgments

J. Samhouri and K. Andrews contributed equally to this work. All authors are especially appreciative of the scientists involved with collection of the data used in these analyses, especially Bill Peterson, the NWFSC trawl survey team, and Sharon Melin. This paper is a result of research supported by the National Oceanic and Atmospheric Administration's Integrated Ecosystem Assessment (NOAA IEA) Program. Contributing support for this work was provided to JS, AS, and MH by the Gordon and Betty Moore Foundation through their generous support of the Ocean Tipping Points project during the development of this paper. This paper is NOAA IEA program contribution no. 2017_2. All co‐authors are grateful for support from the NOAA IEA program and to University of Washington SAFS for providing conference space during the workshop that culminated in this paper, and for socks, particularly the Aquasox, and the Pacific Northwest in general.