High-Resolution Trophic Models Reveal Structure and Function of a Northeast Pacific Ecosystem

Surma, Szymon; Christensen, Villy; Kumar, Rajeev; Ainsworth, Cameron H.; Pitcher, Tony J.

doi:10.3389/fmars.2019.00625

ORIGINAL RESEARCH article

Front. Mar. Sci., 15 October 2019

Sec. Marine Ecosystem Ecology

Volume 6 - 2019 | https://doi.org/10.3389/fmars.2019.00625

High-Resolution Trophic Models Reveal Structure and Function of a Northeast Pacific Ecosystem

$\r\nSzymon Surma*$ Szymon Surma^1*

Villy Christensen¹

Rajeev Kumar²

Cameron H. Ainsworth³

Tony J. Pitcher¹

¹Institute for the Oceans and Fisheries, The University of British Columbia, Vancouver, BC, Canada
²Centre for Fisheries Ecosystems Research, Fisheries and Marine Institute of Memorial University of Newfoundland, St. John’s, NL, Canada
³College of Marine Science, University of South Florida, St. Petersburg, FL, United States

This paper examines the structure and dynamics of the marine ecosystem surrounding Haida Gwaii (an archipelago in the southeastern Gulf of Alaska). Based on previous research, a set of improved mass-balanced models was constructed in Ecopath with Ecosim (EwE) to represent ecosystem states existing circa 1900, 1950, and currently. These models feature increased taxonomic and ecological resolution relative to their predecessors across trophic levels and size classes from zooplankton to whales. A more detailed representation of Pacific herring (Clupea pallasii), including its age structure, predators, and prey, was introduced to permit modeling of the ecosystem role of herring as a forage fish, as well as the ecological impacts of herring fisheries. Gross ecosystem structure and herring trophodynamics were compared across ecosystem states using size spectra and ecological indicators, including mixed trophic impacts (MTIs). The 1950 model was fitted to a comprehensive set of biomass and catch time series. Dynamic ecosystem simulations evaluated the influence of fishing, predation, other natural mortality, and primary productivity trends on ecosystem behavior since 1950, as well as the relative importance of top-down versus bottom-up forcing. Size spectra and ecosystem indicators suggest that the Haida Gwaii ecosystem has not undergone a radical structural shift since 1900, despite heavy exploitation of numerous marine mammals and fish. Moreover, MTIs show that herring constitutes an important mid-trophic level node in the food web, participating in complex interactions with many predators, prey, and competitors. Dynamic ecosystem simulations demonstrate that trends in fishing mortality, trophic interactions, and primary productivity (correlated with the Pacific Decadal Oscillation) are all necessary to explain historical Haida Gwaii ecosystem behavior. These interacting drivers yield a mosaic of top-down and bottom-up trophic control for trophic interactions involving herring and throughout the food web. Simulation results also suggest that production of several herring, salmon, and groundfish stocks may have recently become partially decoupled from primary productivity, perhaps due to changes in copepod guild composition. These results also indicate that a biodiversity decline and “fishing down the food web” occurred off Haida Gwaii since 1950. Finally, the fitted 1950 model provides a robust platform for dynamic ecosystem simulations.

Introduction

This paper aimed to elucidate the structure and dynamics of the marine ecosystem surrounding Haida Gwaii, an archipelago in the southeastern Gulf of Alaska (Figure 1), with particular emphasis on the trophodynamics of Pacific herring (Clupea pallasii). The waters of the study area provide feeding grounds to four herring stocks targeted by commercial and aboriginal fisheries (DFO, 2016). A more detailed description of the oceanography, biogeography, and exploitation history of the productive and diverse Haida Gwaii marine ecosystem can be found in Section A of Supplementary File S1.

FIGURE 1

Figure 1. The study area (∼81,000 km²) located in the southeastern Gulf of Alaska and surrounding the archipelago of Haida Gwaii. Base map obtained from Google Earth Pro. Data LDEO-Columbia, NSF, NOAA. Base map ©2018 Google. Image Landsat/Copernicus. Data SIO, NOAA, U.S. Navy, NGA, GEBCO.

Pacific herring has been classified by Fox et al. (2015, 2018) as a foundation species (sensu Soulé et al., 2003), while Schweigert et al. (2010) consider it the principal forage fish of British Columbia (Canada). In the Northeast Pacific, herring is a notable prey species for numerous commercial fish (Brodeur et al., 2014; Osgood et al., 2016) and seabirds (Schrimpf et al., 2012; Bishop et al., 2015), as well as many marine mammals, including humpback whales Megaptera novaeangliae (Moran et al., 2018; Straley et al., 2018) and Steller sea lions Eumetopias jubatus (Gende and Sigler, 2006; Womble and Sigler, 2006). The trophic role of Pacific herring is especially strong and taxonomically extensive in the spring spawning season, when fish and roe alike are available to a diverse coastal predator guild (Willson and Womble, 2006).

This paper analyzed gross ecosystem structure and Pacific herring trophodynamics (including top-down versus bottom-up control) over time using size spectrum slopes and a suite of ecological indicators. In order to reveal long-term ecosystem shifts, these metrics were compared across three ecosystem states existing circa 1900 (i.e., before industrial whaling), 1950 (i.e., after the main pulse of whaling but prior to that of large-scale commercial herring fishing), and currently. Each ecosystem state was represented by a mass-balanced model built in Ecopath with Ecosim (EwE, Christensen and Walters, 2004), one of the most commonly employed ecosystem modeling frameworks (Coll et al., 2015). The basic principles of EwE are described in Section B of Supplementary File S1. The EwE models presented here, while focusing on Pacific herring, its predators, and prey, exhibit high taxonomic and ecological resolution throughout, with 80 functional groups spanning all trophic levels and size classes from phytoplankton to whales.

This paper also comparatively evaluated the roles of fishing, predation, other natural mortality, and primary productivity as drivers of Haida Gwaii ecosystem dynamics since 1950, as well as the relative importance of top-down versus bottom-up forcing, with particular emphasis on herring trophodynamics. These goals were approached by fitting the 1950 model to an extensive set of biomass and catch time series as well as fishing mortality drivers (1950–2015). The resulting simulated ecosystem dynamics were then statistically related to several oceanographic indices to connect simulation outputs to empirical observations. Furthermore, the model fitting procedure was designed to yield a robust platform for dynamic ecosystem simulations and management strategy evaluation scenarios in the future.

The models employed in this paper derive from a series of previous EwE efforts in the study area. Beattie (1999, 2001) constructed the first Ecopath models of ecosystem states existing off Haida Gwaii (in Hecate Strait) in the early 1900s and 1990s. Vasconcellos and Pitcher (2002) and Ainsworth et al. (2008) expanded the ecological resolution and spatial range of these models to include all northern British Columbia (NBC) waters and developed additional historical models for 1750 and 1950. Ainsworth et al. (2008) also fitted the 1950 NBC model to time series of catch, biomass, and fishing mortality (1950–2000), thus preparing it for use in dynamic ecosystem simulations.

However, the NBC models still featured a relatively simple, aggregated functional group structure at both high (marine mammals, seabirds, elasmobranchs) and low trophic levels (zooplankton, epifaunal invertebrates, and producers). Thus, the detailed investigation of Pacific herring trophodynamics undertaken in this paper necessitated increased group resolution for both non-teleost predators (as generally recommended by Goedegebuure et al., 2017), and zooplanktonic prey of herring. Furthermore, future investigation of the role of herring as a forage fish, and the ecosystem impacts of herring fisheries management, within the research framework outlined by Pikitch et al. (2012) would likewise require higher functional group resolution for herring predators and prey. Recent findings on the great diversity and ontogenetic trajectory of local Pacific herring diet (Pakhomov et al., 2017), representing a great improvement on the simplified herring diet composition in the NBC models, argued for a more detailed multi-stanza representation of juvenile herring age structure. Heterogeneity in biomass and catch trends among local herring stocks (DFO, 2016) suggested the need for their separation in models designed to evaluate ecosystem impacts of herring fisheries. Recent intrusions of warm water masses into the study area, and the likely rise in their frequency with global climate change (Di Lorenzo and Mantua, 2016), prompted the introduction of functional groups for two southern fish species, Pacific hake (Merluccius productus) and Pacific saury (Cololabis saira), which could exert notable pressure on herring through predation and competition, respectively.

In addition, preparation of these models to simulate the likely ecosystem impacts of recolonization of the study area by the sea otter (Enhydra lutris), a keystone species (Estes and Palmisano, 1974) extirpated from Haida Gwaii by the maritime fur trade in the 19th century (Heise et al., 2003), motivated the decision to resolve epifaunal invertebrates and benthic producers in finer ecological detail. Furthermore, the nearly two decades elapsed since 2000 offered the opportunity of adapting the NBC model to a software version (EwE 6) featuring considerable new or improved functionality, including more robust modeling of fish age structure (Christensen et al., 2008). They also allowed the model to be re-fitted to a set of time series expanded and extended with newly available biomass and catch data (e.g., DFO, 2016). Finally, these improvements to the NBC models were designed to avoid the pitfalls of model “recycling” (Essington and Plagányi, 2013). All changes in model structure since the NBC models are summarized in Table 1.

TABLE 1

Table 1. Modifications to functional group structure introduced in the present models.

Materials and Methods

A detailed description of Ecopath model structure and parameterization appears in Section C of Supplementary File S1. Basic Ecopath parameter values for all functional groups and models are given in Section D of Supplementary File S1. Consumer diet compositions appear in Supplementary File S2, while fisheries landings and discards (t km⁻²) by fleet and functional group are given in Supplementary Files S3, S4, respectively. Model parameterization was tested for gross deviations from empirically determined cross-ecosystem patterns using the PREBAL (pre-balancing) technique of Link (2010).

Static Ecosystem Snapshots

Distributions of biomass or abundance versus individual mass or size (commonly known as size spectra) are frequently employed to investigate ecosystem structure and status (Blanchard et al., 2017). When such distributions are plotted using multiple points per size class rather than logarithmically binned classes, and biomasses are not normalized by size interval width, a log–log size spectrum slope of zero is almost universal across marine, freshwater, and terrestrial ecosystems, as revealed by both empirical investigations and models (Blanchard et al., 2017). However, small departures from this value occasionally occur, particularly in heavily exploited ecosystems (e.g., Jennings and Mackinson, 2003).

In order to compare ecosystem structure across all models (1900, 1950, and current), log-transformed functional group biomass densities (t km^–2) were tested for Spearman rank correlation (n = 80) with estimated log-transformed mean individual mass (kg). Where this relationship proved significant (P < 0.05), an analysis of covariance (ANCOVA) was employed to test for differences among model spectrum slopes. Model group biomasses were plotted against individual masses separately, rather than in logarithmically binned size classes, to generate size spectra, and biomasses were not normalized by size interval width. Thus, ANCOVA tested model spectrum slopes for departures from the zero value expected given this approach to defining size spectra (Blanchard et al., 2017). Given the lack of evidence to the contrary (e.g., higher biomasses in 1900 than in 1950 across many size classes), it was hypothesized that size spectrum slopes across all ecosystem states (1900, 1950, and current) would not differ significantly from each other or from zero.

In addition, six ecosystem status indicators were compared across the 1900, 1950, and current models. These included mean trophic levels of the ecosystem (MTL_e) and catch (MTL_c), the Shannon–Wiener biodiversity index (H′) of the ecosystem, total biomass (B_tot) and catch (Y_tot) of non-detritus groups, and total system throughput (TST), i.e., total annual biomass flux (t km^–2 year^–1) through all direct trophic interactions. In the MTL and H′ calculations, functional groups replaced species as the fundamental unit. Based on the findings of Pauly et al. (2001), MTL_e and MTL_c were expected to decline since 1900, as were H′ and TST due to the apparent depletion of many fish and marine mammal groups.

Furthermore, two group-specific indicators were evaluated across the 1900, 1950, and current models. Firstly, proportions of the total herring consumption by predators were calculated for each predator to examine their impacts on herring across ecosystem states. Secondly, mixed trophic impacts (MTIs), measuring the net instantaneous impact (positive or negative; top-down, bottom-up, or competitive) of each functional group on all others (Christensen and Walters, 2004), were examined for all interactions involving herring in the 1950 model. Biomass-weighted average MTI values were computed across all herring stocks and age classes to evaluate overall trophic impacts on and of herring circa 1950. Based on the diverse diet of local herring (Pakhomov et al., 2017) and the importance of this forage fish to numerous predators (Willson and Womble, 2006), MTI values for trophic interactions involving herring were expected to be substantial.

Dynamic Ecosystem Trends

The 1950 Ecopath model was fitted to time series of biomass and catch data (1950–2015) in Ecosim using a non-linear search algorithm manipulating vulnerability parameter values (v), which govern the degree of top-down or bottom-up control of each trophic interaction. Vulnerability values above and below v = 2 imply top-down and bottom-up control, respectively (Christensen et al., 2008). Model fitting followed best practice as outlined by Heymans et al. (2016). The most parsimonious fit to time series among those generated by 10 scenarios was selected using the Akaike Information Criterion with a second-order correction for small sample size (AICc). The 10 fitting scenarios (Table 2) were mainly designed to evaluate the relative importance of trophic interactions and trends in fishing mortality, other natural mortality (M₀), and primary productivity (i.e., bottom-up forcing) as drivers of local ecosystem dynamics through ranking of scenario AICc values. Climatic correlates of the historical primary production anomaly (PPA), a relative trend reconstructed from all time series by the non-linear search algorithm, were investigated using Spearman rank correlations (n = 792) with monthly values of three physical oceanographic indices: the Pacific Decadal Oscillation or PDO (Mantua et al., 1997), the Northern Oscillation Index or NOI (Schwing et al., 2002), and sea surface temperature (SST) at Langara Island off northwestern Haida Gwaii (Figure 1). A full description of the model fitting procedure appears in Section E of Supplementary File S1. Biomass, catch, and fishing mortality time series are provided in Supplementary File S5, while basic Ecosim parameter values appear in Supplementary File S6.

TABLE 2

Table 2. Fitting scenarios tested for the 1950 model.

In order to quantify changing ecosystem structure and status since 1950, trends in three ecological indicators, i.e., mean ecosystem trophic level (MTL_e), the Shannon–Wiener biodiversity index (H′), and the fishing in balance index (FiB), were computed from the most parsimonious fit (i.e., minimum AICc) scenario. FiB is an improved indicator of “fishing down marine food webs” (Pauly et al., 1998) lacking the methodological flaws attributed to MTL_c (Pauly et al., 2000). Significance of observed trends was tested with Spearman rank correlation between time (years) and annual indicator values (n = 66). Based on the findings of Pauly et al. (2001), MTL_e was expected to decline during the fitting interval, while FiB was anticipated to reveal that “fishing down” the local food web had occurred since 1950. A decrease in H′ was also expected based on declines in many functional group biomasses.

In addition, in order to investigate the fullest possible range of trophic and non-trophic factors affecting herring natural mortality M, Spearman rank correlations (n = 26) between herring M and copepod biomass anomalies were computed for 1990–2015. Natural mortalities (DFO, 2016; Figure 2A) for the local (HG, HG 2W, PRD, and CC) and West Coast of Vancouver Island (WCVI) herring stocks were matched with biomass anomalies (Galbraith and Young, 2017; Figure 2B) observed off northwest Vancouver Island for three copepod assemblages (i.e., southern shelf, boreal shelf, and subarctic oceanic). These are distinguished by index species and a descending order of energy content (Mackas et al., 2004). Elevated biomass of energy-rich (i.e., subarctic oceanic and/or boreal shelf) copepods might be expected to decrease starvation mortality in planktivorous fish such as herring, while such mortality could be expected to increase with the biomass of energy-poor (i.e., southern shelf) species (Mackas et al., 2007).

FIGURE 2

Figure 2. Trends in (A) estimated natural mortality (M) for five Pacific herring stocks (DFO, 2016) and (B) biomass anomalies for three copepod assemblages off northwest Vancouver Island, 1990–2015 (Galbraith and Young, 2017).

All statistical tests in the present study were performed using the JMP statistical software package (Altman, 2002). The Bonferroni correction was applied to significance levels for all multiple comparisons.

Results

Static Ecosystem Snapshots

Highly significant (P < 0.001) log–log correlations between functional group biomass and individual mass occurred across all models (1900: Spearman ρ = −0.47; 1950: ρ = −0.54; current: ρ = −0.52). Size spectrum slopes (Figure 3) did not differ significantly among models (ANCOVA: n = 80, df = 2, F = 0.02, P = 0.96), but all departed significantly from the expected value of zero (ANCOVA: n = 80, df = 3, F = 6.69, P < 0.001).

FIGURE 3

Figure 3. Biomass size spectra for Ecopath models representing three ecosystem states (1900, 1950, and current). Functional group biomasses were not normalized by size interval width or grouped into logarithmically binned size classes.

Ecosystem indicator values across all three models appear in Table 3. The mean trophic level of the ecosystem (MTL_e) declined ∼11% between 1900 and 1950 and ∼4% after 1950, while that of the catch (MTL_c) decreased by ∼14 and ∼3% over these intervals. Shannon–Wiener biodiversity (H′) increased ∼4% from 1900 to 1950 and declined ∼3% since 1950. Total living biomass (B_tot) decreased ∼66% between 1900 and 1950 and ∼4% after 1950. Total catch (Y_tot) rose almost sixfold (∼592%) from 1900 to 1950 before falling ∼23% following 1950. Total system throughput (TST) declined ∼80% before 1950 and ∼6% afterward.

TABLE 3

Table 3. Ecosystem indicator values for the three Ecopath models.

The distribution of total herring consumption among predator groups (Figure 4) is a function of proportions of herring in predator diet compositions (Figure 5), predator biomasses, and predator consumption per biomass (Q/B) rates. Herring was generally most prominent in predator diets from 1950, followed by current and 1900 diets (Figure 5A). In 1950, the piscivorous seabird diet included the highest proportion of herring (∼25%), with humpback whale, lingcod (Ophiodon elongatus), and inshore rockfish (Sebastes spp.) diets not far behind (Figure 5B).

FIGURE 4

Figure 4. Proportions of total predator consumption of herring (all stocks and age/size classes combined) attributable to major predator functional groups in the (A) 1900, (B) 1950, and (C) current ecosystem models.

FIGURE 5

Figure 5. Proportions of (A) herring by ecosystem state (1900, 1950, and current) and (B) herring, forage fish (all small pelagic fish except herring, Pacific saury, and eulachon), and euphausiids in the diets of selected predator functional groups.

Large groundfish, primarily Pacific hake and cod (Gadus macrocephalus), arrowtooth flounder (Atheresthes stomias), planktivorous rockfish (Sebastes spp.), and lingcod, were the dominant consumers of herring across all ecosystem states (Figure 4). However, the role of marine mammals was also notable, with cetaceans preeminent among them in 1900 and currently while pinnipeds predominated in 1950. In 1900, humpback and fin whales (Balaenoptera physalus) represented nearly 25% of the total predator consumption of herring, while predatory fish (predominantly Pacific hake, cod, and arrowtooth flounder) accounted for most of the remainder (Figure 4A). In 1950, whale consumption of herring was negligible, with the pinniped (seal and sea lion) share rising to ∼13% of the total and carnivorous fish (mainly Pacific hake, arrowtooth flounder, and planktivorous rockfish) responsible for the rest (Figure 4B). In the current model, pinnipeds and humpback whales accounted for nearly 25% and ∼7% of predator consumption of herring, respectively, while large fish (primarily Pacific hake, arrowtooth flounder, and Pacific dogfish Squalus suckleyi) represented most of the remainder (Figure 4C).

Notable MTIs of herring on other functional groups (MTI > 0.01) are shown in Figure 6A. The number of impacts above this threshold (9) was quite low (∼13%) relative to the number of nodes in the food web (69 if all herring groups are considered together). MTIs of other functional groups on herring are presented in Figure 6B. The number of values exceeding the threshold (21) was higher for impacts on herring (∼30% relative to the number of food web nodes) than that for impacts of herring on other groups. While positive values predominated among the impacts of herring, the opposite was true of impacts on herring. Of the 12 herring predators with MTI > 0.01, 3 registered positive impacts of herring while 9 exerted negative impacts on herring. All three herring prey with MTI > 0.01 showed positive impacts on herring. However, two of these were also negatively impacted by herring, albeit less strongly so. Nine (36%) of the 25 functional groups demonstrating MTI > 0.01 were neither predators nor prey of herring.

FIGURE 6

Figure 6. Notable (>0.01) mixed trophic impacts (MTI) of (A) herring (all stocks and age/size classes combined) on other functional groups and (B) other functional groups on herring in the 1950 ecosystem model. Sm gel stands for small gelatinous.

Dynamic Ecosystem Trends

Akaike information criterion (AICc) values for each 1950 model fitting scenario appear in Table 4. The vulnerability search conducted in the presence of fishing mortality drivers accounted for ∼87% of the AICc reduction in the most parsimonious fit to time series (scenario 10) relative to the baseline (ΔAICc). Fishing mortality drivers alone represented ∼69% of ΔAICc, while a vulnerability search unaided by such drivers achieved only 2–3%. The PPA and non-predation natural mortality forcing functions for selected fish stocks (listed in Section E of Supplementary File S1) accounted for 5 and 3% of ΔAICc, respectively. Running Ecosim with default vulnerability settings (v = 2) for all consumers yielded only ∼3% of ΔAICc.

TABLE 4

Table 4. Adjusted Akaike information criterion (AICc) values and percent improvement in AICc with respect to the scenario 1 baseline (ΔAICc) for each fitting scenario.

The fitted PPA represented a significant increasing trend (i.e., positive correlation with time in years: Spearman ρ = 0.44, P < 0.001) across the fitting interval (1950–2015). It also showed significant positive correlations with local SST (ρ = 0.08, P = 0.02) and the PDO (ρ = 0.26, P < 0.001) and a significant negative correlation with the NOI (ρ = −0.08, P = 0.02). However, correlations with SST and the NOI were no longer significant once the Bonferroni correction for multiple comparisons was applied.

Final vulnerability parameter values appear in Section F of Supplementary File S1. Unlike in the 1950 NBC model fit, no significant correlation was detected between vulnerability and predator trophic level (Spearman ρ = 0.04, P = 0.85). Sixteen consumers, including blue whales (B. musculus), fin, sei (B. borealis), humpback, and sperm whales (Physeter macrocephalus) all adult herring (Figure 7), Pacific cod, and euphausiids, showed v ≥ 2, while 12 others (gray whales Eschrichtius robustus, transient orcas Orcinus orca, seals and sea lions, all seabirds, chinook salmon Oncorhynchus tshawytscha, transient salmon Oncorhynchus spp., eulachon Thaleichthys pacificus, lingcod, and small demersal elasmobranchs) demonstrated v ≤ 2.

FIGURE 7

Figure 7. Biomass time series (squares) and dynamic ecosystem simulation outputs (most parsimonious fit, dotted lines) for adult (age 3+ years) classes of the four local herring stocks.

Significant decreasing trends (i.e., negative correlations with time) were detected for Shannon–Wiener H′ biodiversity (ρ = −0.63, P < 0.001; Figure 8B) and the FiB index (ρ = −0.75, P < 0.001; Figure 8C), but not the mean trophic level of the ecosystem (ρ = −0.11, P = 0.30; Figure 8A). The Shannon–Wiener H′ showed a particularly steep decline after 2005. The FiB index oscillated between positive and negative values early in the fitting period but became consistently negative by 1970.

FIGURE 8

Figure 8. Trends in the (A) ecosystem mean trophic level MTL_e (non-significant), (B) Shannon–Wiener biodiversity H′ (significant), and (C) fishing in balance index FiB (significant) from dynamic ecosystem simulation outputs (most parsimonious fit).

The statistical analysis (Table 5) of trends (Figure 2) in copepod assemblage biomass anomalies (Galbraith and Young, 2017) and herring natural mortalities (DFO, 2016) yielded the expected negative correlation between subarctic oceanic copepod biomass and herring M for the CC and WCVI stocks. However, while this correlation was significant at α = 0.05, it became marginally non-significant for CC with the Bonferroni correction for multiple comparisons. Similar correlations for HG and HG 2W M were already marginally non-significant at α = 0.05. Also as expected, southern shelf copepod biomass exhibited a positive correlation with HG M (α = 0.05). However, this correlation too became marginally non-significant under the Bonferroni correction. More surprisingly, boreal shelf copepod biomass showed positive correlations with M (PRD, CC, and WCVI) that remained significant under the Bonferroni correction.

TABLE 5

Table 5. Spearman rank correlations between natural mortalities for five British Columbia herring stocks (Haida Gwaii, Haida Gwaii 2 West, Prince Rupert District, Central Coast, and West Coast of Vancouver Island) and abundance anomalies for three copepod assemblages (southern shelf, boreal shelf, and subarctic oceanic) off northwestern Vancouver Island (1990–2015).

Discussion

Static Ecosystem Snapshots

The size spectrum slope of approximately −0.20 seen across all models (Figure 3) differed subtly but significantly from the expected value of zero predominant in almost all ecosystem categories (Blanchard et al., 2017). However, the slope observed in this study was almost identical to that detected in the North Sea by Jennings and Mackinson (2003), although the latter only sampled size classes between zooplankton and medium-sized fish. These deviations from the expected slope of zero may have derived from predator–prey mass ratios increasing with individual predator mass (e.g., baleen whales versus zooplankton) or from underestimated abundance of the highest size classes (i.e., whales and sharks). The lack of significant differences in size spectrum slope among models supported the hypothesis that biomass distributions across size classes were broadly similar in all three ecosystem states (1900, 1950, and current). It also indicated that the deviations from the expected spectrum slope of zero observed in all models are unlikely to stem from heavier exploitation of large size classes, as such exploitation only began in earnest after 1900 with industrial whaling (Webb, 1988) and pinniped hunting (Olesiuk, 2010, 2018). In the late 19th century, at least several large fish groups (e.g., halibut, salmon, and dogfish) had also yet to be heavily exploited in the study area (Chittenden, 1884). Furthermore, the biomasses of several groups in small to mid-size classes, from phytoplankton to small fish and squid, had to be increased in the 1900 model to satisfy the direct and indirect trophic demands of large, unexploited whale populations. This steepened the biomass size spectrum slope, and suggested the possibility of higher ecosystem productivity before commercial whaling than indicated by the NBC models (Ainsworth et al., 2008).

The modest reductions in MTL_e and MTL_c across the three ecosystem states (Table 3) also did not suggest drastic shifts in trophic structure since 1900. The reductions in these two indices likely resulted from depletion of upper-level predators (i.e., marine mammals and large groundfish) during the 20th century. The unexpected peak in H′ observed for 1950 may have reflected declines in several mid-trophic level groups (e.g., small squid, forage fish, and herring) after 1900. High biomasses and P/B or Z-values for such mid-trophic groups, as well as zooplankton and phytoplankton were necessary, as mentioned above, to satisfy direct and indirect demands of predation by unexploited whale populations in 1900. Depletion of whales (Surma and Pitcher, 2015), pinnipeds (Olesiuk, 2010, 2018), and large predatory fish, and the resulting relaxation of their trophic demands, apparently underlay the strong declines in B_tot and TST between 1900 and 1950. The increase in Y_tot between 1900 and 1950 presumably reflected greater fishing power due to the advent of internal combustion engines and improved (e.g., powered) gear, while the decrease in Y_tot after 1950 resulted mainly from the closure of the massive herring reduction fishery in 1963 following the depletion of all local herring stocks (Hourston, 1978).

The distribution of herring consumption among predators and across ecosystem states (Figure 4) resulted from population trends in humpback and fin whales (abundant in 1900, depleted by 1950, currently recovering; Surma and Pitcher, 2015) and pinnipeds (depleted in 1900 and 1950, currently abundant; Olesiuk, 2010, 2018). Schweigert et al. (2010) observed a pattern in groundfish versus marine mammal consumption of herring off western Vancouver Island similar to that detected in the present study. However, neither the composition of total herring consumption by predators (Figure 4) nor MTIs (Figure 6B) truly supported the hypothesis of top-down impacts of Pacific cod on herring (Walters et al., 1987) or the synchronous changes in consumption of herring by all predators suggested by Oken et al. (2018) based on an older Hecate Strait EwE model.

The notable role of herring in 1950 predator diets relative to the 1900 and current models (Figure 5A) resulted from a higher ratio of herring to other prey biomass in 1950, yielded by the historical diet calculation algorithm described in Section C of Supplementary File S1 (Buszowski and Steenbeek, 2017). While herring comprised 15–25% of several marine mammal, seabird, and large fish diets in 1950, no evidence of predator specialization on herring was apparent (Figure 5B). However, this may have been an artifact of functional group aggregation, with group diet compositions, e.g., biomass-weighted averages of constituent species diets for seabirds, obscuring species-specific dietary requirements.

The MTIs of herring (Figure 6A) and on herring (Figure 6B), combined with the distribution among predators of the total consumption of herring (Figure 4) and the proportions of herring in predator diets (Figure 5), indicated that interactions between herring and its predators are frequently asymmetric (i.e., more strongly impact either predator or prey). Inshore rockfish are an exception to this rule, featuring nearly equal MTI of herring (bottom-up, Figure 6A) and on herring (top-down, Figure 6B). Similarly, two herring prey groups (macrozooplankton and amphipods) exerted large positive MTI on herring (Figure 6B) but experienced only slightly lower negative MTI (Figure 6A). The intricate combination of top-down and bottom-up interactions apparently connecting herring with its predators and prey (Figure 6) suggested that the “wasp-waist” trophic control mechanism proposed for many forage fish, in which the latter impact predators and prey through bottom-up and top-down control, respectively (Cury et al., 2000), most likely underestimates the complexity of herring trophodynamics in the modeled food web. Herring also shared mutually negative MTI with forage fish, indicating competition, while transient orcas indirectly benefited herring through cascading effects of predation on smaller marine mammals (i.e., small odontocetes and pinnipeds). Phytoplankton exerted a strong positive impact on herring, suggesting indirect bottom-up control via zooplankton. These static Ecopath results assisted in framing expectations for simulations of herring trophodynamics in Ecosim.

Dynamic Ecosystem Trends

The ecosystem simulations revealed that while fishing mortality, trophic interactions (governed by vulnerability values), and primary productivity were all necessary to emulate ecosystem dynamics since 1950 and fit the model to time series, none of these factors alone was sufficient (Table 4). These findings accorded with those derived from several previous EwE models of Northeast Pacific ecosystems, including the Gulf of Alaska (Gaichas et al., 2011), southeast Alaska, and the Aleutian Islands (Heymans et al., 2005). However, as in the 1950 NBC model (Ainsworth et al., 2008), fishing mortality trends exerted an especially strong influence on the most parsimonious fit to time series (and hence on ecosystem dynamics), while the effect of primary productivity was much weaker (Table 4). This pattern has been observed in a diverse, global array of modeled marine ecosystems in both temperate and tropical waters, although its opposite occurred in several ecosystems (Mackinson et al., 2009). Furthermore, a spectrum stretches from fully fishing-dominated to climate-dominated marine ecosystems (Gaichas et al., 2011), and model fits improved by primary productivity (though not necessarily dominated by it) are even more ubiquitous than those influenced by fishing mortality (Mackinson et al., 2009). In addition, the reduction in AICc due to the reconstructed primary productivity trend in the present 1950 model (5%) was over twice that (2%) in its NBC counterpart.

The primary production anomaly (PPA) reconstructed from all available time series was not particularly similar to the input phytoplankton biomass series, as shown in the last panel of Supplementary Figure S1.3 in Supplementary File S1. This contrast was also observed in an EwE model of the English Channel ecosystem (Araújo et al., 2006). However, the increasing trend in the PPA after 2000 agreed with satellite measurements of surface chlorophyll a density in the study area (O’Brien et al., 2017). Intriguingly, correlations between the PPA and all three physical oceanographic indices (i.e., local SST, PDO, and NOI) in the present 1950 model showed reversed signs relative to those observed in the corresponding NBC model (Ainsworth et al., 2008). The significant positive correlation between the PPA and PDO in the present model agreed with that detected in a southeast Alaska ecosystem model fitted to time series from 1963 to 2001 (Heymans et al., 2005), although Guénette (2005) obtained no such correlation for a similar interval. The correlations found in the present study and by Heymans et al. (2005) may be explained by the “optimal stability window” hypothesis (Gargett, 1997), according to which positive PDO values increase water column stability, a moderate level of which supports optimal primary productivity in the Northeast Pacific. The dome-shaped relationship between water column stability and primary productivity (Gargett, 1997) may explain the weakness of the observed correlation between PPA and PDO. However, this correlation remains consistent with the “optimal stability window” hypothesis, since large positive monthly PDO values (presumably indicating periods of strong water column stability) are quite rare (Mantua et al., 1997). If “optimal stability windows” correspond to moderately positive PDO values, the weak positive correlation between PPA and PDO observed in this study is as could be expected. Local and regional SST was positively correlated with the PPA and PDO (Mantua et al., 1997), respectively, although the correlation between the local SST and PPA lost its significance under the Bonferroni correction for multiple comparisons. Nevertheless, the observed correlations between the PPA, PDO, and SST supported the role of SST as a driver of fisheries production in the Gulf of Alaska detected by Mcowen et al. (2014).

However, while the PPA assisted the model in emulating overall ecosystem dynamics since 1950, it also yielded unobserved and apparently spurious recoveries of several fish functional groups (all herring, coho salmon O. kisutch, chinook salmon, Pacific cod, lingcod, and sablefish Anoplopoma fimbria), adversely affecting their fits to time series. These were improved, along with AICc (Table 4), by forcing functions for non-predation natural mortality (Supplementary Figure S1.6 in Supplementary File S1). This finding suggested that production of several local fish stocks might have become at least partly decoupled from local primary productivity since 2000.

Climate-associated fluctuations in the relative biomass of three copepod assemblages (Figure 2B) differing in energy content (subarctic oceanic, boreal shelf, and southern shelf) suggested a potential mechanism for such decoupling. Results of the correlation analysis described above (Table 5) indicated that while natural mortality trends in several herring stocks (Figure 2A) might be associated with shifts in the relative biomasses of copepod assemblages (and thus with average copepod energy content), these associations remain to be conclusively demonstrated. Low sample size (n = 26) due to the short time series available (1990–2015) most likely reduced the statistical power of the correlation analyses. In addition, many extraneous factors, including herring stock assessment model structure and smoothing of natural mortality time series derived from these models (DFO, 2016), could have confounded these putative correlations. Possible mesoscale variability in copepod assemblage biomass anomalies between the waters northwest of Vancouver Island (where the copepod data were obtained) and the study area (where the modeled herring stocks feed) might also have affected correlations for these stocks. This inference was supported by the expected negative correlation between herring natural mortality and the subarctic oceanic copepod biomass anomaly showing unambiguous significance only for the WCVI stock, whose feeding range overlaps the copepod sampling area.

It is, however, noteworthy that hypothesized correlations between southern shelf and subarctic oceanic copepod biomass anomalies, respectively, and herring natural mortalities were clearest for the two modeled herring stocks most strongly depressed in recent decades (i.e., HG and CC). While the significantly positive correlations between the boreal copepod biomass anomaly and many herring stock (PRD, CC, and WCVI) natural mortalities were somewhat surprising, several boreal shelf copepod species were rare in, or absent from, recent local herring diets (Pakhomov et al., 2017). These positive correlations could also derive from trophic effects not captured in the present ecosystem models (e.g., elevated boreal copepod biomass increasing the productivity of herring competitors or predators).

The incomplete dataset on copepod assemblage biomass anomalies in Hecate Strait showed similar interannual and decadal patterns as its full counterpart for northwestern Vancouver Island (Galbraith and Young, 2017), as did abundance data for large versus small copepods from the central Gulf of Alaska (Batten et al., 2018). This indicated that shared dynamics likely affect copepod community composition (and thus average energy content) in much of the Northeast Pacific. This inference was supported by significantly positive correlations linking both abundance and northward extent of southern copepods to the PDO (Batten and Walne, 2011). These, in turn, helped explain the apparent synchrony between high herring natural mortalities (Figure 2A) and positive PDO values between late 2002 and late 2007. Copepod biomass anomalies were also associated with Northeast Pacific coho salmon marine survival and growth rates, as well as sablefish recruitment (Mackas et al., 2007). Thus, these anomalies could explain the discrepancies between observed and expected Ecosim biomass trends for coho and sablefish, which were largely eliminated by forcing non-predation natural mortality (Supplementary Figure S1.6 in Supplementary File S1). For coho and chinook salmon, competition with pink salmon (O. gorbuscha) might also partially explain the observed lack of recovery (Ruggerone and Irvine, 2018).

The fact that the net contribution of the PPA to the most parsimonious fit to time series exceeded that of the trophic control (vulnerability) parameters (Table 4) indicated a substantial role for bottom-up forcing in Haida Gwaii ecosystem dynamics since 1950. This inference was echoed by the low vulnerability values obtained for 12 functional groups, suggesting bottom-up control of interactions between these groups and their prey (Christensen and Walters, 2004; Christensen et al., 2008). These findings agreed with the overall bottom-up control of Northeast Pacific fish production detected by Ware and Thomson (2005). However, the high vulnerability values obtained for 16 functional groups (including all adult herring and most large cetaceans, e.g., humpback whales) indicated top-down control. These results largely accorded with MTIs of and on herring observed in the static 1950 model (Figure 6). They also supported the hypothesized top-down control of prey populations by whales in the North Pacific (Croll et al., 2006), particularly humpback whales preying on herring in the Gulf of Alaska (Heintz et al., 2010; Surma and Pitcher, 2015; Moran et al., 2018; Straley et al., 2018).

The significant decrease in Shannon–Wiener H′ biodiversity (Figure 8B) over the fitting interval (1950–2015) likely reflected increased evenness in biomass distribution among functional groups. This, in turn, resulted from depletion of several fish stocks abundant in 1950, as well as partial recovery of initially depleted whale populations (Supplementary Figures S1.2, S1.1 in Supplementary File S1, respectively). The significant decline in the FiB to consistently negative values since 1970 (Figure 8C) indicated that “fishing down the food web” had occurred in the study area. This finding agreed with the conclusion of Pauly et al. (2001) based on MTL_c trends for Canadian Pacific fisheries catches but conflicted with that of Fu et al. (2012), which relied on an older EwE model of the Hecate Strait ecosystem. However, the decrease in FiB observed in the present study was not accompanied by a significant decline in MTL_e (Figure 8A), likely due to simultaneous marine mammal (i.e., whale and pinniped) recoveries and herring declines (Figure 7). The MTL_e trend detected in the present analysis was relatively similar to that observed by Fu et al. (2012), albeit shifted downward by approximately one trophic level.

Fits to biomass time series obtained in this study for the major herring stocks (HG, PRD, and CC; Figure 7) and Pacific cod (Supplementary Figure S1.2 in Supplementary File S1) were similar to those generated in an OSMOSE model of the NBC food web (Fu et al., 2017). This suggested the potential for multimodel inference, perhaps within a formal framework such as that developed by Spence et al. (2018). It was also noteworthy that herring fits to time series were improved by using fishing mortality drivers and catch time series based on ratios of reconstructed to reported catch from the Sea Around Us project database (Ainsworth, 2018). This suggested that true catches for HG, HG 2W, PRD, and CC herring in the 1950s and 1960s likely exceeded values utilized in stock assessments.

Caveats and Limitations

It must be admitted here that while every attempt was made to employ the most complete and reliable data available and follow best practice in the model parameterization and fitting processes, the latter were necessarily burdened with many forms of error. Firstly, the increased functional group resolution of the present models could have raised total parameter uncertainty. However, the static and dynamic model parameter values were constrained within plausible limits by the pre-balancing tests designed by Link (2010) and by fitting to an expanded and extended set of time series, respectively, both of which strongly mitigated the potential increase in total uncertainty. The pre-balancing checks confirmed that static parameter values in the present models do not depart strongly from those expected based on empirically determined cross-ecosystem patterns. This suggested that parameter uncertainty lies within acceptable limits. In addition, the impacts of static parameter uncertainty were limited by the fact that in dynamic ecosystem simulations, these values served merely as starting points from which model runs were initialized. Furthermore, the augmented functional group resolution of the present models reduced the structural uncertainty generated by aggregating multiple ecologically diverse species (e.g., seabirds) within a single group with common parameter values. Indeed, even in the present models, the multispecies composition of numerous functional groups could have masked specialist predator responses to biomass dynamics of individual prey species (e.g., herring), ultimately affecting the performance of dynamic ecosystem simulations and projections of overall ecosystem behavior (Pinnegar et al., 2005). In addition, the present models still feature rather simplified representations of lower trophic levels, which could have masked the complex impacts of bottom-up forcing, including oceanographic regime shifts (Mantua et al., 1997; Schwing et al., 2002), on ecosystem dynamics.

Secondly, the absence of explicit age structure and stock-recruitment responses for most piscivorous fish in the present models could have obscured the effects of strong trophic ontogeny in such predators on their interactions with forage fish (Essington et al., 2015; Walters et al., 2016; Hilborn et al., 2017). Thirdly, the non-linear search algorithm fitting vulnerabilities and primary production anomalies to time series can generate satisfactory fits to data with model outputs not representative of the actual ecosystem, since multiple output configurations can yield similarly parsimonious fits (Ainsworth et al., 2008). The search algorithm may also become trapped in local minima well above the true minimum AICc. Fourthly, structural and parameter error in assessment models affect their outputs, which provide static parameter values as well as biomass and mortality time series.

Fifthly, fear-mediated, sublethal effects (Frid et al., 2008; Wirsing et al., 2008) of top predator (e.g., Pacific sleeper shark, Somniosus pacificus) presence on mesopredator (e.g., harbor seal, Phoca vitulina) foraging behavior, in which predator avoidance by mesopredators yields cascading effects on prey (e.g., Pacific herring), are absent from the present models, although their ecosystem impacts could be substantial. Sixthly, mass-balanced model outputs could noticeably underestimate the importance of energy-rich forage fish (e.g., herring and eulachon) to many predators relative to energy-balanced models (Surma et al., 2018). Finally, non-trophic factors not represented in these food web models (e.g., chemical and microplastic pollution, parasitism and infectious disease, ocean warming, acidification, and hypoxia) may also noticeably affect food web function, necessitating a cautious evaluation of the importance of fishing mortality and primary productivity as drivers of ecosystem dynamics (Mackinson et al., 2009). Nevertheless, pre-balancing tests and fitting the 1950 model to time series using the non-linear search algorithm greatly reduced the likelihood of erroneous representations of ecosystem dynamics, explanations of past trends, and projected responses to future perturbations.

Conclusion

The present models provide an improved representation of the ecosystem states existing off Haida Gwaii and in the southeastern Gulf of Alaska in 1900, 1950, and currently. While the Ecopath model structure focuses on Pacific herring, its predators, and prey, it also features increased functional group resolution relative to its predecessors across all trophic levels and size classes from zooplankton to whales. This improved resolution facilitated the investigation of interactions between herring and its predators and prey in greater detail than was possible previously, though both static ecosystem analyses (e.g., MTIs) and dynamic ecosystem simulations examining the ecological role of herring since 1950. The results of this investigation suggested that herring likely plays an important role as a mid-trophic level node (forage fish) in the Haida Gwaii food web by participating in a complex network of top-down and bottom-up interactions with numerous predators and prey as well as competing with other small planktivorous fish.

Comparative analysis of the three present models indicated that ecosystem structure (as reflected in biomass size spectrum slope) did not undergo radical change since 1900, despite heavy exploitation of many marine mammal and fish groups. The 1950 model, fitted to time series following best practices, revealed that fishing mortality, trophic interactions, and primary productivity correlated with oceanographic indices are all required to explain historical dynamics in the Haida Gwaii ecosystem, generating a combination of top-down and bottom-up trophic control. It also indicated that the production of some herring, salmon, and groundfish stocks may have recently been partially decoupled from primary productivity, perhaps due to oceanographically driven changes in copepod guild composition. Dynamic ecosystem simulations further showed that a biodiversity decline and “fishing down the food web” occurred off Haida Gwaii between 1950 and 2015. Finally, the fitted 1950 model provides a robust platform for future dynamic ecosystem simulations and management strategy evaluation scenarios, particularly those designed to investigate the ecosystem impacts of herring fisheries.

Data Availability Statement

All datasets generated for this study are included in the manuscript/Supplementary Files.

Author Contributions

SS was responsible for parameterization of most new functional groups, all historical models, biomass spectra, model fitting, statistical analyses, and writing of the manuscript. VC designed the model fitting procedure and suggested introducing the whale biomass accumulation rates. RK ported the latest NBC model to EwE 6, as well as parameterizing the hake group and herring stock and age structure in the current model (in the historical models this was done by SS). CA developed the original NBC models and prepared some of the time series used in model fitting. TP supervised the work of SS, RK, and CA, and oversaw the general direction of model development. All authors edited the manuscript.

Funding

SS received a doctoral scholarship (PGS-D) partially supporting this work from the Natural Science and Engineering Research Council of Canada (grant number 6564). TP received a Strategic Grant partially supporting this work from the Natural Science and Engineering Research Council of Canada (grant number 447247). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Conflict of Interest

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.

Acknowledgments

We thank Evgeny A. Pakhomov (Director, Institute for the Oceans and Fisheries, UBC and Professor, Department of Earth, Ocean, and Atmospheric Sciences, UBC) and Andrew Trites (Director, Marine Mammal Research Unit, Institute for the Oceans and Fisheries, UBC) for contributing novel herring diet data and expert opinion on local whale abundances, respectively; and colleagues from the Pacific Biological Station (Fisheries and Oceans Canada) for assistance in obtaining necessary data. We specifically thank Jaclyn Cleary, R. Ian Perry, and John K. B. Ford for providing herring stock assessment model outputs, zooplankton biomass anomaly values, and transient orca abundance data, respectively. We also thank Jeroen Steenbeek and Joe Buszowski (EwE Research and Development Consortium) for their repeated technical assistance, including preparation of custom EwE functionality; Dr. Divya Varkey (Fisheries and Oceans Canada) and Dr. Mimi Lam (Centre for the Study of the Sciences and the Humanities, University of Bergen) for supporting model development; Dr. Wade Smith (Institute for the Oceans and Fisheries, UBC), Jeff Scott (M.Sc., Institute for the Oceans and Fisheries, UBC), Sahir Advani (Ph.D. candidate, Institute for the Oceans and Fisheries, UBC), and Nicolás Talloni-Alvarez (Ph.D. candidate, Institute for the Oceans and Fisheries, UBC) for their general encouragement and support; the Haida people and the non-Haida inhabitants of Haida Gwaii for their interest in this research; and specifically thank Russ Jones (Nang Jingwas), hereditary chief of the Naa Saagaas XaaydaGaay clan and retired head of the Haida Oceans Technical Team. Finally, we acknowledge the prior presentation of this research in the Ph.D. dissertation of the corresponding author (Surma, 2019).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2019.00625/full#supplementary-material

References

Ainsworth, C. (2018). British Columbia marine fisheries catch reconstruction, 1873 to 2011. BC Stud. 188, 81–89. doi: 10.14288/bcs.v0i188.186115