This study examined four Large Marine Ecosystems (LMEs: Alaska-Eastern Bering Sea, California Current, Northeast US and northern Gulf of Mexico) that are part of NOAA's IEA program (Levin et al., 2009, 2014). Each LME is a distinct type of marine ecosystem (Figure 1) and these LMEs collectively represent not only latitudinal, bathymetric, productivity, and exploitation gradients, but also encapsulate a wide range of variable habitats and taxa groups useful for contrasts (Murawski et al., 2010). The Eastern Bering Sea (hereafter Alaska), is a sub-arctic, high productivity system that is characterized by an extensive gradually sloping shelf and a deep sea basin. While the system has a relatively low population density compared to other IEA LMEs in the US, it experiences a high level of human activity in terms of commercial fishing, and to a lesser extent oil and gas development and transport (Zador et al., 2014). The California Current is a temperate Pacific coast ecosystem that extends from southern British Columbia, Canada to Baja California, Mexico. The California Current LME produces abundant ecosystem goods and services including fisheries, recreation, tourism, and energy production. This ecosystem is fueled by seasonal upwelling of cold nutrient water resulting in a very productive system (Levin et al., 2013). The Northeast US is a temperate Atlantic ecosystem with high productivity that supports a diverse array of invertebrates, pelagic fish, groundfish, seabirds, and marine mammals. With a history of fisheries spanning centuries, the region has experienced sustained impacts on the marine ecosystem (Link et al., 2002). The northern Gulf of Mexico (hereafter Gulf of Mexico) is a semi-enclosed, sub-tropical coastal sea in the Atlantic that supports a large recreational and commercial fishing industry and also provides many goods and services such as, oil and gas production, tourism, and habitat for endangered species (Karnauskas et al., 2013). Further descriptions of these ecosystems are provided in the Supplementary Materials.

The total model prediction performance from the gradient forest analysis (the proportion of variance explained in a random forest) averaged across the suite of ecological indicators from each LME ranged from 0.01 to 0.07 (R²_s; Table 3). The model prediction performance for the cumulative ecological indicators had a range of 0–0.30 (R²_s), which is consistent with other studies using similar analyses by Large et al. (2015b; 0–0.21 R²_s), Pitcher et al. (2011; 0–0.35 R²_s), but lower than Baker and Hollowed (2014; 0–0.77 R²_s). These model prediction values may appear low; however, in studies with similar results using aggregated ecological indicator data (containing more unexplained variability than raw species abundance metrics) the gradient forest was still able to identify important variables and thresholds (Large et al., 2015b). In contrast to these past studies, this particular analysis focused on cross-ecosystem comparisons and only included ecological indicators that were transferrable across the differing LMEs into the analysis. Inclusion of LME-specific ecological indicators (e.g., longhorn sculpin in the Northeast US; Methratta and Link, 2006) improved the variance explained, but would not allow for true cross-ecosystem comparisons, thus they were excluded from this particular study and account for the lower overall variance explained. Additionally, and more to the point in this comparative study, although the variance explained may be low, multiple (but very different) statistical analyses yielded similar threshold points both within and across ecosystems. In all ecosystems, the ecological indicators that were included in the gradient forest models were mean length (structural) and richness and diversity (resilience) indicators (Table 3). In the Atlantic ecosystems the models maintained a higher number of the ecological indicators in the gradient forest analysis with the fewest in the Alaska ecosystem.

Both the gradient forest analyses and GAMs did not identify a single driver that was consistently dominant across the four ecosystems, though fisheries landings was an important predictor in models for all systems (Figure 2, Table 4). Generally, anthropogenic variables ranked high in their impact, explaining a larger portion of the changes in the ecosystem trends (at least three of the five most important pressure variables were anthropogenic for all four regions in the gradient forest analysis) than environmental pressures. GAM relationships between ecosystem trends and pressure variables also showed that a larger number of anthropogenic variables significantly impacted ecosystems compared to environmental variables, given the relatively higher values of deviance explained (Table 4). Trends of pressure variable importance were consistent with the historical understanding of each region, with the Northeast US and California Current being strongly impacted by the anthropogenic pressures assessed. The pressures that were specific to a region (e.g., ice cover, hypoxic area; Table 2) were ranked less important than the large scale climatic pressures such as, PDO and AMO in explaining patterns of the cumulative ecological indicators. Exploitation had a high impact in the Atlantic regions.

Common landings and exploitation thresholds of ~2–4.5 t km⁻² landings and ~20–40% exploitation of the total estimated biomass were detected with the gradient forest analyses (Figures 3A,C), bearing in mind that the relative importance of exploitation was low in Alaska and the Gulf of Mexico (Figure 2). Individual thresholds of ecological indicators showed how structural, functional and biodiversity properties of each ecosystem shifted across individual pressure gradients (Figure 4; see Supplement 2). Lagged fishing pressures (landings and exploitation) showed very similar threshold mean and ranges to their non-lagged counterparts (Figures 3B,D). This was also confirmed by the GAM of ecosystem trends (Table 4). The thresholds identified by the GAMs of ecosystem trends along pressure gradients were similar to the thresholds identified by the gradient forest analysis, indicating these ecosystem-level thresholds are robust (Figure 5; see Supplementary Materials). Ecosystem thresholds in response to yearly human population increases were lowest in coastal communities that were (already) most densely populated (Northeast US, Gulf of Mexico); in the lowest populated regions (Alaska) the mean thresholds were highest and they had the widest range (Figure 3E). Yearly changes in the GDP showed ecosystem shifts occurring at higher GDP increases in the Northeast US, Gulf of Mexico, and California Current, compared to Alaska where large ranges of GDP were observed (Figure 3F). SST showed differing thresholds between regions that were unique to the specific ecosystems, largely reflective of the latitudinal position of each ecosystem (Figure 3G). Thresholds for ecological indicators occurred in more positive phases of the AMO, PDO, and MEI (Figures 3H,I), although these tended to have a wide range of threshold. Freshwater anomalies generally ranked lower in terms of importance in explaining ecosystem shifts, but showed a wide threshold region in the California Current (Figure 3J). Chlorophyll a had a higher and wider threshold range in California Current (Figure 3K), indicating that higher levels of primary production in terms of chlorophyll concentration is an important driver there, likely due to upwelling. The three other regions indicate a basal concentration of ~0.7 mg m⁻³ chlorophyll a to avoid ecosystem shifts (Figure 3K).

The gradient forest analysis identified the relative size of the cold pool in Alaska as the only region specific pressure variable that ranked within the top five important pressures in explaining ecosystem shifts (Figure 2). The ecosystem thresholds occurred when the cold pool was ~0.2–0.3x relative to the previous year. Ecosystem thresholds along ice retreat index gradients were ~10–50%, while thresholds along north Pacific index were ~−2.0 to 0. In the California Current, ecosystem thresholds were around ~−1.0 to 1.0 for mean annual north Pacific gyre oscillation, 3,000 m³ s⁻¹ 100 m⁻¹ total upwelling magnitude index and northern oscillation index at ~−4. Region specific thresholds in the Northeast US were between ~−2.0 and 2.0 for winter north Atlantic oscillation, ~ 0.06–0.07 N m⁻² for wind stress and ~1.0–1.3 GS. Ecological indicator thresholds along gradients of mean current transport in the Gulf of Mexico had the lowest importance for the regions according to the gradient forest analysis (Figure 2). Hypoxic area (threshold at ~1.5–2.0 × 10⁴ km²) and Atlantic warm pool (threshold at ~2.0–2.5 × 10⁶ km²) in the Gulf of Mexico explained 0.01 and 0.012 of the ecological indicator response, respectively; which was relatively low compared to other pressure variables in that LME.

Our results demonstrate that there are consistent patterns in ecosystem response from common pressures across four large marine ecosystems, and despite multiple potential mechanisms, the detected trends and thresholds to such pressures in these ecosystems were remarkably repeatable. Although each ecosystem examined has different socio-economic histories (Hollowed et al., 2011; Link et al., 2011a; Karnauskas et al., 2013; Levin et al., 2013), different levels of population density, and differential reliance on living marine resources that vary in the use of marine habitats and ultimately shape the stability of the ecosystems, there are a number of common trends that are surprisingly consistent across all ecosystems. These common trends would be difficult to detect if not examined at the ecosystem-level and in a multivariate context, both in terms of detecting baseline reference points and observing emergent properties of marine ecosystems. That these common trends and thresholds exist is insightful for both further understanding of marine ecosystems, as well as management thereof.

One key result is that, at an ecosystem-level, removals of biomass (via landings-based exploitation) do have repeatable and consistent thresholds. There are different ecological mechanisms in which such ecosystem-level responses can be observed, but consistently there is an impact to overall size, congruent with overfishing theory, as well as tendencies toward smaller organisms with hyper-exploitation (Pauly et al., 1998; Pauly and Palomares, 2005; Shackell et al., 2010; Darimont et al., 2015; Worm and Paine, 2016). This exploitation impacts the biomass composition in an entire ecosystem, shifting either biodiversity or measures of biomass ratios (e.g., pelagic to demersal) and implies that exploitation selectively impacts certain facets of an ecosystem consistent with what is known for the ecological effects of targeted fisheries (Shin et al., 2005; Jennings and Collingridge, 2015). Additionally, some form of broad-scale climate forcing is regularly identified as a key driver of ecosystem-level responses. The ecological mechanisms for this can vary, but largely filter through changes in vital rates and related bioenergetics processes (Holsman et al., 2016).

Another commonality was that when examining the impacts of the pressure variables to cumulative ecosystem responses, anthropogenic pressures rank high. Due to the pressure-response relationship of these analyses, this does not necessarily reflect the current status of a given ecosystem, but rather implies that certain pressures have heavily impacted these ecosystems within the history of the time series analyzed. This is not to say that environmental pressures are not important; rather that the anthropogenic pressures tended to more consistently emerge as clearer features that impact observed ecosystem dynamics.

As a particular example of such anthropenic pressures, in all the ecosystems, landings have decreased since the 1970s and 1980s, and some fish stocks have experienced a phase of rebuilding (Rosenberg et al., 2007; Worm et al., 2009; Lotze et al., 2011). Ecosystem responses to commercial fishing are surprisingly consistent in the four regions examined in this study. Ecosystem thresholds were observed at landings of ~2–4.5 t km⁻² and fisheries exploitation of ~20–40% of the total estimated biomass. All landings estimates fall within previously determined ecosystem-level surplus production model threshold estimates of 1–6 t km⁻² (Bundy et al., 2012; Link et al., 2012; Lucey et al., 2012), further corroborating the robustness of a potentially universal ecosystem-level fisheries yield, at least in the northern hemisphere.

Lower population increases and higher GDP in coastal communities were related to ecosystem shifts where population density is highest (population density of Northeast US: ~300 indv km⁻²). These threshold values increased with decreasing population densities (Gulf of Mexico: ~56 indv km⁻², California Current: ~50 indv km⁻². Alaska: ~0.5 indv km⁻²; U.S.Census Bureau, 2016). Human migration toward more heavily populated areas appears to have a disproportionately large impact on ecosystems. This could relate to a longer history of commercial fishing in a particular region or issues with urban infrastructure (e.g., sewage treatment, erosion prevention). Likely, pressures relating to human population increase act cumulatively to account for these patterns in ecosystem shifts (Halpern et al., 2008; Stallings, 2009; Madin et al., 2016). The Northeast US, California Current, and Gulf of Mexico regions exhibited higher annual GDP increases which were related to ecosystem shifts, compared to Alaska where ecosystem shifts occurred when GDP was both increasing and decreasing. There was a wider range where ecosystem shifts could occur along gradients of population and GDP increase in Alaska, suggesting that less infrastructure, fewer types of industry or lower climate regulation influence ecosystem dynamics and stability (Li and Fang, 2014).

In terms of environmental pressures, another commonality was evidence that all ecosystems appear to influenced by multi-annually varying climate drivers, albeit seen via different indices in each region. The discovery of large-scale climate patterns have been an important step in connecting climate to ecological, biological, and oceanographic patterns in marine ecosystems (Mantua and Hare, 2002; Alheit et al., 2014). Of the environmental drivers examined in this study, the large-scale climate drivers (PDO and AMO) were important in explaining ecosystem shifts when examined at a yearly time scale. Seasonal patterns of climate pressures, however, were not examined in this study to maintain consistency with other indicator and pressure variables, but other relationships and thresholds may certainly exist at different time scales. Both positive PDO and AMO anomalies are associated with dry, hot temperatures in the U.S. (McCabe et al., 2004) and generally correlate positively with SST anomalies that subsequently can cause a shallower ocean mixed layer and lower primary productivity (Mantua and Hare, 2002; Nye et al., 2014). When PDO and AMO anomalies appear to be positive, the ecosystem shifts toward a threshold response. Patterns of chlorophyll across all four ecosystems indicate a base level of primary productivity needed for sustainable fisheries at chlorophyll a concentrations of ~0.7 mg m⁻³, noting that chlorophyll concentration is positively associated with fisheries yields (Friedland et al., 2012). In the California Current, ecosystem shifts were observed at higher concentrations of chlorophyll which are linked periods of strong upwelling, but ecosystem stability is maintained at lower chlorophyll concentrations (Kahru et al., 2012). Scatterplots of monthly fisheries yields vs. chlorophyll concentration by Friedland et al. (2012) show spikes in observed fisheries yields when chlorophyll concentrations were at or above ~0.7 mg m⁻³.

The analyses in this study were not used to examine specific mechanistic links between specific pressures and responses, but rather to identify and compare significant threshold ranges of ecosystems (represented by indicators) along pressure gradients. Individual ecological indicator thresholds can also be examined across individual pressure gradients to determine specific reference points at the single indicator level (e.g., Samhouri et al., 2010; Large et al., 2013). This allows for the specific examination of shifts across gradients of individual ecological indicators that can then be examined using other methods (such as, Generalized Linear Models) to determine a trend or directional shift. While doing so would further elucidate the dynamics of a given ecosystem, this study is primarily aimed at examining the location of ecosystem thresholds across pressure gradients across ecosystems. We acknowledge that further insight into individual indicator trends and thresholds are important to gain a fuller understanding of a given ecosystem. We also assert that comparisons across ecosystems, that then detect common trends and thresholds, are equally powerful.

Mechanistic links have been made between fisheries production and multiple drivers including fishing, trophodynamics, and the environment (Gaichas et al., 2012a,b; Holsman et al., 2012; Link et al., 2012; Pranovi et al., 2014; Longo et al., 2015). Multivariate thresholds analyses enable scientists and policy makers to present the complexity of multiple pressures on whole ecosystems in a way that is easier to communicate and comprehend (Peterman, 2004; Large et al., 2015a,b). The methods presented here are complex, in terms of statistical analyses, but are able to offer clear ecosystem-level outputs of threshold ranges across pressure gradients. As such, thresholds provide a powerful tool to delineate and communicate quantifiable tipping points for ecosystems (Foley et al., 2015).

Each ecosystem in this study has, at some stage in the last half century, experienced overfishing (Bakkala et al., 1979; Rosenberg et al., 2007; Walters et al., 2008; Link et al., 2011a; Miller et al., 2014). The impacts of overfishing are complex, but well researched. For example, it has been found that the concentration of fishing on large, predatory species can destabilize ecosystems, resulting in a loss of biodiversity and increasing blooms of lower-trophic organisms, often undesirable species like jellyfish (Pauly et al., 1998; Purcell et al., 2007). More recently, the position of humans within food-webs has been categorized as “hyperkeystone” (Worm and Paine, 2016) or as “super predators” (Milius, 2015), indicating that the current impacts of humans on food webs (marine or otherwise) could lead not only to lower abundances of predatory species, but also size-stunted predator populations with an impaired ability to regulate prey species even in rebuilding scenarios (Darimont et al., 2015). Here we present ecosystem thresholds that can be used as reference points to support coordinated efforts to develop sustainable fisheries via EBM policies that support rebuilding strategies for depleted fish stocks, ecosystem aggregate yield limits, and to explore the social tradeoffs and potential social benefits of changing how we use living marine resources (Murawski, 2000; Balmford, 2002; Howarth and Farber, 2002; Rosenberg et al., 2007; Worm et al., 2009; Khan and Neis, 2010; Link, 2010; Plagányi et al., 2014; DePiper et al., in press).

The thresholds developed here can also be used to build proactive strategies to avoid regime shifts due to overfishing, population increase and climate change, particularly when explored through simulation modeling (Samhouri et al., 2010; Fulton et al., 2011; Fay et al., 2013, 2015; Large et al., 2013). In particular, given the strength of bottom up controls on the systems we evaluated, climate change may be of increasing impact. There have been major breakthroughs in understanding how long-term climatic events and climate change will impact marine species. Many marine communities are predicted to have continued impact by changes in mean SST and are predicted to respond in numerous ways from range shifts (Pinsky and Fogarty, 2012; Gattuso et al., 2015; Heenan et al., 2015) to phenological shifts (Parmesan, 2006; Brown et al., 2016). While the outcomes of climate change are predicted to be overwhelmingly negative, there may be some economic and social benefits for well-managed fisheries in warming scenarios (Barange et al., 2014). Warming climate events in the Gulf of Alaska during the 1970s and 1980s caused dramatic shifts in catch compositions of groundfish, and actually increased the catch biomass of higher trophic-level groundfish by over 250% (Anderson and Piatt, 1999). Anchovy and sardines are known to have multidecadal regime shifts in the Pacific Ocean due to natural, long-term climate variability (Chavez et al., 2003). This can drastically change the catch compositions of these ecosystems over long periods of time and it is important to understand the consequences of unintentionally overfishing during these temporary regime shifts for the long-term sustainability of marine ecosystems. Systems approaches to management, including EBM, can contribute to successful resilience of ecosystems by improving the ability to detect and react to ecological feedbacks (Hughes et al., 2005). Climate thresholds presented in this study can be incorporated into local resource and ocean management and used to develop rules around total harvests of specific species based on known or predicted distributional shifts (Link et al., 2011b; Pinsky et al., 2013; Heenan et al., 2015).

While the location of ranges of ecosystem-level thresholds along both human and environmental pressure gradients are easily interpreted, these insights are best made against a backdrop of dynamic biological and environmental conditions. While it is easy for many to agree that benefits to human well-being correlate positively with ecosystem services, it is often difficult to incorporate these ideas into management (Arkema et al., 2015). Social indicators that examine the mechanistic links between human well-being and ecosystems are being examined, but are generally underdeveloped (McShane et al., 2011; Colburn and Jepson, 2012; Howe et al., 2014; Pollnac et al., 2015; Yang et al., 2015). Studies such as the one undertaken here address issues from an ecosystem-based fisheries management standpoint, but the inclusion of more social indicators (such as, GDP and human population increases) would move this work more fully toward EBM. Improving and expanding on indicators for human well-being would greatly enhance the ability to examine thresholds for societal needs within the context of this study, but also to examine the inherent tradeoffs between the needs of people and marine ecological sustainability and conservation (Dearing et al., 2014; DePiper et al., in press).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.