De ﬁ ning ecosystem thresholds for human activities and environmental pressures in the California Current

. 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 suf ﬁ ciently to pro-voke major responses of species, habitats, and communities. We developed a framework based on multimodel inference to de ﬁ ne ecosystem-based thresholds for human and environmental pressures in the California Current marine ecosystem. To demonstrate how to apply the framework, we explored two dec-ades 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 identi ﬁ ed the existence of threshold responses of ﬁ ve 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 modi ﬁ cation, and (2) sea lion pup production and the summer mode of the Paci ﬁ c Decadal Oscillation (PDO). Considered collectively, however, these alternative analytical approaches imply that as many as ﬁ ve of the nine ecosystem states may exhibit threshold changes in response to negative PDO values in the summer (copepods, scavengers, ground ﬁ sh, and marine mammals). This result is consistent with the idea that the in ﬂ uence of the PDO extends across multiple trophic levels, but extends current knowledge by de ﬁ ning 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 miti-gation may enhance progress toward ecosystem-based management goals.

INTRODUCTION Any performance management systemwhether 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 1999). 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. 2004, Arkema et al. 2006, Levin et al. 2009, Samhouri et al. 2013, Link and Browman 2014, a frequent first step has been to identify and report on indicators of ocean conditions (Kershner et al. 2011, Andrews et al. 2015. In marine EBM and elsewhere, however, indicators are necessary but not sufficient for successful performance management (Otley 1999). Useful indicators must describe a system in relation to the objectives set for it (Melnyk et al. 2004). 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 , Hamel et al. 2015, UMCES 2016. 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 1992). One of the most informative techniques to establish reference points relies upon thresholds in ecosystem state. Following Groffman et al. (2006), 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 1998, Methratta and Link 2006, Martin et al. 2009.
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. 2008). 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 2006, Litzow and Hunsicker 2016. Many analytical techniques for identifying ecosystem thresholds based on state-pressure relationships have also been proposed (Samhouri et al. 2010, Large et al. 2013, 2015a, b, Hunsicker et al. 2016, 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.oceaneco nomics.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. 2014, Leising et al. 2014, ONMS 2015, Zador 2015. 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 ❖ www.esajournals.org 2 June 2017 ❖ Volume 8(6) ❖ Article e01860 in many ecosystems ) and can alter the human-derived benefits of a system . As such, we believe our MMI framework that identifies ecosystem thresholds has broad application to other regions in the process of operationalizing EBM.

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. 2011, Scheffer et al. 2012). The workflow (Fig. 1) can be summarized in four parts, including (1) pre-treatment, (2) screening, (3) functional form identification, and (4) threshold identification.
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 preestablished indicator frameworks (e.g., European  Zuur et al. 2003). 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. 2014). 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. 2010, Large et al. 2013, 2015a, b, Baker and Hollowed 2014. 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. 2010, Large et al. 2013, 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. 2009, Foley et al. 2015.

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. 2014). 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. 2009. 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. 2009, Andrews et al. 2015. 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. 2013) and the Northern Oscillation Index (NOI; Schwing et al. 2002). 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. 2008), and decadal changes in sea surface temperature regimes in the northeast Pacific are tracked by the Pacific Decadal Oscillation (PDO; Mantua and Hare 2002). 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. 2011), including shifts in planktonic communities  and lower trophic-level fishes (Chavez et al. 2003), as well as higher trophiclevel fishes (Lindley et al. 2009), seabirds, and mammals (Leising et al. 2014).
Human activities within the CCS are also quite dynamic (Andrews et al. 2015). Economic shocks (e.g., McKenna et al. 2012), emerging technologies (Kim et al. 2012, Plummer andFeist 2016), and regulatory shifts (e.g., in fisheries; Hilborn et al. 2012, Lubchenco et al. 2016) 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. 2009). Here, we drew from the 19 human activities presented in the CCIEA ) and focused on the 10 that provided available data across most of the period of interest (Table 1).
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 Note: Further details of data sets can be found in Appendix S1: Table S2.
❖ www.esajournals.org 5 June 2017 ❖ Volume 8(6) ❖ Article e01860 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 2; Williams et al. 2014). 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. 2014 for details).
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. 2003, Garc ıa-Reyes and Largier 2012), we focused on six basin-scale environmental pressures in the North Pacific (Table 2; Appendix S1: Table S2). These environmental pressures are the primary basin-wide indicators of oceanography and climate in the CCIEA ) 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. 2011).
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. 2008 and regional (Halpern et al. 2009, Andrews et al. 2015 attention in relation to changes in marine ecological integrity (Table 1; 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. 2015 for details).  Williams et al. (2014). 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. ❖ www.esajournals.org 6 June 2017 ❖ Volume 8(6) ❖ Article e01860 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 2002, Ellis et al. 2012, Baker and Hollowed 2014, Large et al. 2015b. 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. 2014), especially given the limitations of inferences based on P-values (Barber and Ogle 2014).
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 2), but the span of overlap among all ecosystem and pressure time series was 10 yr (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(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 2001, Large et al. 2015b. 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 2 -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 2 -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 2 -importance that exceeded 0.01.
For pressures that increased cumulative ecosystem response R 2 -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 2 -importance values of each ecosystem state given a split at specific values of the pressure (Large et al. 2015b). We used the R packages "randomForest" (Liaw and Wiener 2002, R Core Team 2016) and "gradientForest" (Ellis et al. 2012) 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. 2009, Bestelmeyer et al. 2011). 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. (2012) 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 ❖ www.esajournals.org 7 June 2017 ❖ Volume 8(6) ❖ Article e01860 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. 1-3).
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. 1, 2 and 4). Formally, let E y be the value of the ecological indicator in year y, a a fixed intercept effect, P y the value of the pressure in year y, s() the smoothing function, e 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: The formulation of Eq. 3 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 q. In contrast, the formulation of Eq. 4 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. 2009); (2) the generalized cross validation score (GCV; Wood 2004b) 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 2002). For both GAM and GAMM, we used thin plate regression spline smoothing terms, the "mgcv" (Wood 2004a) package, and set the size of the basis dimension to 4 to reduce the possibility of overfitting (Large et al. 2013).
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. 2000, Bestelmeyer et al. 2011, Large et al. 2013. 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 ❖ www.esajournals.org 8 June 2017 ❖ Volume 8(6) ❖ Article e01860 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. 1 (except for pre-treatment steps described in Methods). 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 pressurestate relationships based on a single method (Figs. 2-4, Table 3; Appendix S1: Table S3). 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. 2, Table 3; 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 3; 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. 3 and 4). 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. 3). 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. 4).
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. 4). 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. 3). 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. 4). The contrast in the functional form of this pressure-state relationship demonstrates the potential for strong ❖ www.esajournals.org 9 June 2017 ❖ Volume 8(6) ❖ Article e01860 SAMHOURI ET AL.

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 3, Fig. 5; 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).
For the thresholds related to human pressures (Table 3, Fig. 5; 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. 4, Table 3). At the high end, the winter northern copepod anomaly changed by 180% across the threshold found with NPGO w (Fig. 4, Table 3).

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. 6). 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. 6a). 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. 6b).

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€ om et al. 2009, Steffen et al. 2011, Scheffer et al. 2012. In marine environments, the consequences of crossing thresholds-for population extinctions, anoxia, and acidification, for instance-are often pressures using data from 2003 to 2012. Top: importance of environmental and human pressures weighted across all ecosystem states. Pressures with R 2 -importance ≥0.01 were considered capable of predicting variation in ecosystem states. Cumulative importance (in R 2 -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. (Fig. 2. Continued) ❖ www.esajournals.org 11 June 2017 ❖ Volume 8(6) ❖ Article e01860 difficult to reverse (deYoung et al. 2008. 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. 2006). 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. 2011). 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. 5). 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. 4). 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 Fig. 3. Truncated GAM analyses (data from 2003 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.
❖ www.esajournals.org 12 June 2017 ❖ Volume 8(6) ❖ Article e01860 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 3). 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. 2014), 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. 2015a), though we did identify threshold changes in the cumulative ecosystem response using gradient forest analysis (Fig. 2). 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. 2015). 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. 2015).  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 3).

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. 2008, Bograd et al. 2009, Black et al. 2011. Previous studies have established the strong influence of human activities in this system (Halpern et al. 2009, Hilborn et al. 2012, Andrews et al. 2015. 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. 2017).
Nonetheless, the analyses presented here stand in contrast to results from other systems. For instance, in the Northwest Atlantic fisheriesrelated pressures were associated with substantial threshold responses in the ecosystem (Large et al. 2013(Large et al. , 2015b. 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. 4). 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. 2008), 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. 4). 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, Fig. 6. 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. 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 2003). While the importance of PDO in the California Current has a long history (Hare et al. 1999, Mantua andHare 2002), evidence for associated nonlinear ecosystem responses is more limited. For example, Bestelmeyer et al. (2011) 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 2016). 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. 5a, Table 3), a result consistent with the idea that the influence of the PDO extends across multiple trophic levels (Hare et al. 1999, Mantua andHare 2002). The location of thresholds associated with PDO w for copepods and marine mammals ranged more widely (Fig. 5b, Table 3), 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. 2).
Recognition of taxon-specific vs. ecosystemwide 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. 6). 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. 2) that spanned the same range of values as thresholds for individual ecosystem states identified via  and gradient forest analysis (Figs. 2 and 5). 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. 2009), 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 ❖ www.esajournals.org 16 June 2017 ❖ Volume 8(6) ❖ Article e01860 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. 2013(Large et al. , 2015bTam 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.   shaded area is last five years of dataset, arrow shows whether the indicator has increased (↗), 24 decreased (↘), or remained the same (→) over the last five years, symbol shows whether the 25 mean of the last five years of the dataset was > (+), < (-) or within (•) 1 s.d. of the long-term 26 mean. See Table 2 for data sources. 27 28 Figure S2. Time series of basin-scale environmental pressure indicators for the California 29

LITERATURE CITED
Current ecosystem from 1996 -2014. Symbology same as Figure S1. 30 Figure S3. Time series of coastwide human pressure indicators for the California Current 32 ecosystem from 1996 -2014. Symbology same as Figure S1. 33 34 Figure S4. Truncated GAM analyses (data from 2003 showing ecosystem response (top 36 panels), 1st derivative (s′(X); middle panels) and 2nd derivative (s″(X); bottom panels) to 37 environmental or human pressures, where the dashed black line is the GAM smoother, gray 38 polygon is 95% CI, points are raw data, thick solid line indicates the threshold range where the 39 95% CI of the 2nd derivative does not include 0 (bold polygons on 2nd derivative plots), and red 40 dotted arrow indicates the best estimate of the location of the threshold (i.e., where the 2nd 41 derivative is at its absolute maximum value within the threshold range). (Note that each top panel  42 is identical to a panel in Figure 3). i) Figure S5. Full GAM analyses (data from 1996 showing ecosystem response (panels 48 a,d,g,j), 1st derivative (s′(X); panels b,e,h,k) and 2nd derivative (s″(X); panels c,f,i,l) ecosystem 49 responses to environmental or human pressures, where the dashed black line is the GAM 50 smoother, gray polygon is 95% CI, points are raw data, thick solid line indicates the threshold 51 range where the 95% CI of the 2nd derivative does not include 0 (filled polygons on 2nd 52 derivative plots), and red dotted arrow indicates the best estimate of the location of the threshold 53 (i.e., where the 2nd derivative is at its absolute maximum value within the threshold range). 54 (Note that panels labeled (a,d,g,f) are identical to panels in Figure 4).