# Spatial and temporal dynamics of predator-prey species interactions off western Canada

Caihong Fu, Norm Olsen, Nathan Taylor, Arnaud Gruss, Sonia Batten, Huizhu Liu, Philippe Verley, Yunne-Jai Shin

<figure><img src="image_1.png" /></figure>

## To cite this version:

Caihong Fu, Norm Olsen, Nathan Taylor, Arnaud Gruss, Sonia Batten, et al., Spatial and temporal dynamics of predator-prey species interactions off western Canada. ICES Journal of Marine Science, 2017, 74 (8), pp.2107 - 2119. 10.1093/icesjms/fsx056. hal-01928592

HAL Id: hal-01928592

https://hal.umontpellier.fr/hal-01928592

Submitted on 8 Jun 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ICES Journal of Marine Science (2017), 74(8), 2107–2119. doi:10.1093/icesjms/fsx056

Original Article

Spatial and temporal dynamics of predator-prey species
interactions off western Canada

Caihong Fu,<fnref n="1" />* Norm Olsen,<fnref n="1" /> Nathan Taylor,<fnref n="1" /> Arnaud Grüss,<fnref n="2" />,<fnref n="3" /> Sonia Batten,<fnref n="4" /> Huizhu Liu,<fnref n="5" /> Philippe Verley,<fnref n="6" /> and Yunne-Jai Shin<fnref n="6" />,<fnref n="7" />

Fu, C., Olsen, N., Taylor, N., Grüss, A., Batten, S., Liu, H., Verley, P., and Shin, Y.-J. Spatial and temporal dynamics of predator-prey species interactions off western Canada. – ICES Journal of Marine Science, 74: 2107–2119.

Received 12 January 2017; revised 9 March 2017; accepted 15 March 2017; advance access publication 12 May 2017.

Ecosystem models are valuable tools for informing fisheries management due to their ability to simulate the spatial dynamics of modelled species, their trophic interactions, and their responses to fishing in an ecosystem context. In this study, we developed an OSMOSE (Object-oriented Simulator of Marine Ecosystems Exploitation) model for the Pacific North Coast Integrated Management Area (PNCIMA) ecosystem off western Canada, which simulated the entire life cycle of six key species and for the first time integrated spatial population structure and “background” taxa. Background taxa are of secondary importance for the study at hand but have the potential to be important prey or predators of the key species. The primary aim of the study was to explore how the population dynamics of the key species differed over time and different management areas, with results focusing on Pacific Herring (*Clupea pallasii*) and Pacific Cod (*Gadus macrocephalus*) that have been assessed on a single-species basis in the last 5 years. Results found that the population dynamics of a specific species varied in different management areas due to differences in species interactions particularly in the form of predation mortality, which supports the current area-specific assessment and management framework. The study also indicated that increasing predation mortality may have caused the decline of the Pacific Cod populations. By contrast, increasing starvation mortality was found to be a limiting factor for the Pacific Herring populations. The discoveries from these OSMOSE simulations provide important information for fisheries management within the PNCIMA ecosystem.

**Keywords:** ecosystem modelling, fisheries management, OSMOSE, predation mortality, spatial population structure, species interaction.

Introduction

Quantitative ecosystems models, i.e. mathematical representa-
tions of ecosystem dynamics, are valuable tools for studying

complex interactions and their effects (Seidl, 2017). Species inter-
actions and their spatial patterns are fundamental to the dynam-
ics of fish populations and ecosystems, though often ignored in

<footnote n="1">¹Fisheries and Oceans Canada, Pacific Biological Station, 3190 Hammond Bay Road, Nanaimo, BC V9T 6N7, Canada</footnote>

<footnote n="2">²Department of Marine Biology and Ecology, Rosenstiel School of Marine and Atmospheric Science, University of Miami, 4600 Rickenbacker Causeway, Miami, FL 33149, USA</footnote>

<footnote n="3">³Sustainable Fisheries Division, Southeast Fisheries Science Center, 75 Virginia Beach Drive, Miami, FL 33149-1099, USA</footnote>

<footnote n="4">⁴Sir Alister Hardy Foundation for Ocean Science, c/o 4737 Vista View Cr, Nanaimo, BC V9V 1N8, Canada</footnote>

<footnote n="5">⁵Department of Computing Science, Vancouver Island University, Nanaimo, BC V9T 5S5, Canada</footnote>

<footnote n="6">⁶Institut de Recherche pour le Développement (IRD), UMR MARBEC 248, Centre de Recherche Halieutique Méditerranéenne et Tropicale, Avenue Jean Monnet, CS 30171, 34203 Sète Cedex, France and Université de Montpellier, Place Eugène Bataillon, CC093,34095 Montpellier cedex 5, Bâtiment 24, 34095, France</footnote>

<footnote n="7">⁷Marine Research (MA-RE) Institute and Department of Biological Sciences, University of Cape Town, Private Bag X3, Rondebosch 7701, South Africa</footnote>

<footnote>*Corresponding author: tel: +1 250 7298373; fax: +1 250 7567053; e-mail: caihong.fu@dfo-mpo.gc.ca.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Map of the PNCIMA showing the spatial grid of the PNCIMA OSMOSE model (light grey cells), the three management areas (delineated by double dashed lines) of Herring PRD, HG, and CC, and the two management areas (separated by dashed lines) of Cod and Lingcod (Hecate Strait and Queen Charlotte Sound).</figcaption></figure>

single-species-based fisheries management. Travis *et al.* (2014) argued that fisheries science and management must develop a sharper focus on species interactions, and understand the consequences of disruption of these interactions. Individual-based models (IBMs) that integrate driving ecological processes at individual, population, community, and ecosystem levels, are powerful tools for predicting ecosystem dynamics in a changing world that can provide valuable information to managers (Grimm *et al.*, 2017). OSMOSE (Object-oriented Simulator of Marine Ecosystems Exploitation) is such a spatially explicit IBM that accounts for both size-based trophic interactions and whole-life-cycle dynamics of marine organisms (Shin and Cury, 2004). OSMOSE has been employed to model trophic dynamics and the impacts of fishing in a variety of ecosystems, for example, the Southern Benguela (Shin and Cury, 2004; Travers-Trolet *et al.*, 2014), the West Florida Shelf (Grüss *et al.*, 2016a), and the Strait of Georgia in Canada (Fu *et al.*, 2013). However, spatial differences in species interactions and their implication in fisheries management have never been explicitly investigated with OSMOSE.

Because simulations with OSMOSE necessitate extensive information on entire life cycles, typically no more than 10 to 15 key species are included. By limiting the number of species included in an OSMOSE model, the computation time and memory capacity can be kept reasonable while the complex interactions within the study ecosystem are simplified. However, this limited number of species can cause the OSMOSE model to miss important, and even sometimes major, prey and predators when the study ecosystem is characterized by a relatively high biodiversity (e.g. the West Florida Shelf, Grüss *et al.*, 2016a). In the present study, we enhanced OSMOSE to allow explicit consideration of nearly all the taxa of a given ecosystem without compromising much of the

computation time and without requiring extensive information on whole life cycles. This enhancement was done through the inclusion of “background” taxa, that is, taxa that are of secondary importance for the study but have the potential to be important prey or predators of the modelled key species.

We applied the enhanced OSMOSE model to the ecosystem of the Pacific North Coast Integrated Management Area (PNCIMA) off western Canada to investigate species interactions and their spatial pattern in relation to population spatial structure. We addressed two fundamental questions related to species interactions, specifically in the form of predation mortality. (i) What were the levels of predation mortality for the different spatial populations of the modelled key species? (ii) Which species and taxa contributed to the predation mortality and how did they differ for the different spatial populations of the same key species?

## Material and methods

### Study area and OSMOSE model components

The PNCIMA off western Canada encompasses ~88 000 km², extending from the coastal watersheds to the outer limit of the continental slope (Figure 1). The PNCIMA is bounded to the north by the Canada–Alaska border and to the south by Brooks Peninsula on northwest Vancouver Island and Quadra Island to the east of Vancouver Island (Lucas *et al.*, 2007). The PNCIMA ecosystem supports hundreds of species, a few dozen of which are commercially exploited (Lucas *et al.*, 2007). The spatial grid of the OSMOSE model developed in the present study covers the whole PNCIMA region, divided into grid cells of 10 × 10 km² (Figure 1).

We focused on six key species in the PNCIMA OSMOSE model: Pacific Herring (*Clupea pallasii*), Pacific Cod (*Gadus*
<table><caption>**Table 1.** Key species, background taxa and LTL groups included in the PNCIMA OSMOSE model.</caption><thead><tr><th>Species/Taxon name</th><th>Type</th><th>Species represented</th></tr></thead><tbody><tr><td>Euphausiids</td><td>Key species</td><td>Thysanoessa spp., Euphausia spp.</td></tr><tr><td>Pacific Herring</td><td>Key species</td><td>Pacific Herring (Clupea pallasi)</td></tr><tr><td>Arrowtooth Flounder</td><td>Key species</td><td>Arrowtooth Flounder (A. stomias)</td></tr><tr><td>Walleye Pollock</td><td>Key species</td><td>Walleye Pollock (T. chalcogramma)</td></tr><tr><td>Pacific Cod</td><td>Key species</td><td>Pacific Cod (G. macrocephalus)</td></tr><tr><td>Lingcod</td><td>Key species</td><td>Lingcod (O. elongatus)</td></tr><tr><td>Seals and sea lions</td><td>Background taxa</td><td>Harbour Seal (Phoca vitulina), Steller Sea Lion (Eumetopias jubatus)</td></tr><tr><td>Whales</td><td>Background taxa</td><td>Humpback Whale (M. novaeangliae)</td></tr><tr><td>Pacific Hake</td><td>Background taxa</td><td>Pacific Hake (M. productus)</td></tr><tr><td>Pacific Ocean Perch</td><td>Background taxa</td><td>Pacific Ocean Perch (Sebastes alutus)</td></tr><tr><td>Spiny Dogfish</td><td>Background taxa</td><td>Spiny Dogfish (S. acanthias)</td></tr><tr><td>Flatfish</td><td>Background taxa</td><td>Dover Sole (Microstomus pacificus), Rock Sole (Lepidosetta bilineata), English Sole (Parophryys vetulus), Sand Sole (Psettichthys melanostictus), Rex Sole (Glyptocephalus zachirus), Flathead Sole (Hippoglossoides elassodon), Butter Sole (Isopsetta isolepis), Curlfin Sole (Pleuronichthys decurrens), and Starry Flounder (Platichthys stellatus), Pacific Sanddab (Citharichthys Sordidus), Yellowfin Sole (Limanda Aspera), Slender Sole (Lyopsetta Exilis), C-O Sole (Pleuronichthys Coenosus)</td></tr><tr><td>Pacific Halibut</td><td>Background taxa</td><td>Pacific Halibut (H. stenolepis)</td></tr><tr><td>Petrale Sole</td><td>Background taxa</td><td>Petrale Sole (E. jordani)</td></tr><tr><td>Shelf Rockfish</td><td>Background taxa</td><td>Yellowtail (Sebastes flavidus), Silvergray (Sebastes brevispinis), Widow (Sebastes entomelas), Bocaccio (Sebastes paucispinis), Canary (Sebastes pinniger)</td></tr><tr><td>Slope Rockfish</td><td>Background taxa</td><td>Yellowmouth (Sebastes reedi), Roughey (Sebastes aleutioanus), Redstripe (Sebastes proriger), Sharpchin (Sebastes zacentrus), Redbanded (Sebastes babcocki), Shortspine Thornyhead (Sebastolobus altivelis), Splitnose (Sebastes diploproa), Longspine Thornyhead (Sebastolobus alascanus), Darkblotched (Sebastes cremeri), and Shortraker (Sebastes borealis)</td></tr><tr><td>Inshore Rockfish</td><td>Background taxa</td><td>Yelloweye (Sebastes rubberrimus), Quillback (Sebastes maliger), Copper (Sebastes caurinus), China (Sebastes nebulosus), Tiger (Sebastes nigrocinctus)</td></tr><tr><td>Spotted Ratfish</td><td>Background taxa</td><td>Spotted Ratfish (Hydrolagus colliei)</td></tr><tr><td>Sablefish</td><td>Background taxa</td><td>Sablefish (Anoplopoma fimbria)</td></tr><tr><td>Coho Chinook</td><td>Background taxa</td><td>Coho Salmon (Oncorhynchus kisutch), Chinook Salmon (Oncorhynchus tshawytscha)</td></tr><tr><td>Shallow Benthic Fish</td><td>Background taxa</td><td>Eelpouts (Zoarcidae), Poachers (Agonidae), Sculpins (Cottidae)</td></tr><tr><td>Forage Fish</td><td>Background taxa</td><td>Eulachon (Thaleichthys pacificus), Smelts (Osmeridae), Sandlance (Ammodytes hexapterus)</td></tr><tr><td>Crabs</td><td>Background taxa</td><td>Dungeness crab (Metacarcinus magister), Tanner Crab (Chionoecetes spp.), Red Rock Crab (Cancer productus)</td></tr><tr><td>Shrimp</td><td>Background taxa</td><td>Smooth Shrimp (Pandalus jordani), Spiny Shrimp (Pandalus borealis eous), Pink Shrimp (Pandalus goniturus), Sidestrip Shrimp (Pandalopsis disbar), Prawn (Pandalus platycterus)</td></tr><tr><td>Detritus Benthos</td><td>Background taxa</td><td></td></tr><tr><td>Phytoplankton</td><td>LTL group</td><td></td></tr><tr><td>Copepods</td><td>LTL group</td><td></td></tr></tbody></table>

The reference species of each background taxon (in bold) is used for obtaining mean length and weight at different stages as well as distribution maps.

macrocephalus), Lingcod (Ophiodon elongatus), Arrowtooth Flounder (Atheresthes stomias), Walleye Pollock (Theragra chalcogramma), and Euphausiids (Thysanoessa spp. and Euphausia spp.). These six key species are typically treated as resident species and have potentially strong trophic interactions with one another and the other species of the PNCIMA ecosystem (Pearsall and Fargo, 2007). Within the PNCIMA, Pacific Herring (hereafter abbreviated as Herring) has been assessed and managed as three separate populations in three distinct areas: Prince Rupert District (PRD), Haida Gwaii (HG), and Central Coast (CC) (Figure 1). Genetic microsatellite variation (Beacham et al., 2008) and long-term tagging studies (Flostrand et al., 2009) identified these Herring populations as distinct, likely due to differences in timing of spawning as well as geographic isolation of spawning populations. Similarly, Pacific Cod (hereafter abbreviated as Cod) and Lingcod have been assessed and managed as separate populations in two distinct areas: Hecate Strait (HS) and Queen Charlotte Sound (QCS) (Figure 1). We therefore chose to treat Herring, Cod, and Lingcod within each management area as individual populations. Subsequently, a total of 10 populations were

individually modelled in OSMOSE (i.e. Herring-PRD, Herring-HG, Herring-CC, Cod-HS, Cod-QCS, Lingcod-HS, Lingcod-QCS, Walleye Pollock, Arrowtooth Flounder, and Euphausiids). The PNCIMA OSMOSE model also includes 19 background taxa and 2 lower-trophic-level (LTL) groups (Table 1). We used the functional groups in Ainsworth (2006) as background taxa unless local experts (e.g. G. Workman, pers. comm.) disagreed, as was the case with rockfish. The rockfish species were combined according to their spatial distribution rather than their feeding guild (Ainsworth, 2006). The novel inclusion of background taxa allowed us to explicitly consider all species and taxa that constitute potential predators and prey of the modelled key species. A diagrammatic representation of the hierarchical architecture of key species, background taxa and LTL groups in the PNCIMA OSMOSE model and of the interactions between these three categories of marine organisms is provided in Figure 2.

PNCIMA OSMOSE simulates the life cycles of the 10 populations, from the egg stage to the terminal age, at a time step of three months. At the first time step following the production of eggs, the total number of eggs of each population is split into 120
super-individuals called “schools” (Figure 2), which are distributed spatially according to input distribution maps. The distribution maps (10 × 10 km²) are density-based and obtained from geo-referenced data of both commercial fisheries and research surveys (data archives being maintained by Fisheries and Oceans Canada, Pacific Biological Station, Nanaimo, British Columbia). At each time step, OSMOSE simulates the biological and ecological processes of these schools, including growth, predation, starvation, “diverse” natural mortality ($M^{diverse}$) due to causes unaccounted for by the model, fishing, reproduction, and spatial movement (including migration).

The average growth of these schools follows the von Bertalanffy growth model, but the growth rate of a specific school is

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Diagrammatic representation of the different categories of marine organisms considered in the PNCIMA OSMOSE model and of the interactions between these categories. The three categories of marine organisms considered in PNCIMA OSMOSE include: key species, background taxa, and LTL groups. The direction of the arrows represents trophic flow. The dynamics of key species are simulated on an individual basis, while that of background taxa are simulated on a biomass basis. LTL groups simply serve as prey for key species and background taxa. Key species are divided into several populations, and populations are further divided into schools, which consist of individuals that have the same age, body size, food requirements and, at a given time step, the same spatial coordinates. The total biomass of each background taxon is apportioned into three life stages (young-of-the-year, juveniles, and adults), converted to abundance, and then apportioned to 120 schools.</figcaption></figure>

determined by its consumption rate and prey availability (Shin and Cury, 2004). Predation is assumed to be an opportunistic process, occurring under the conditions of size suitability (within a minimum and a maximum predator to prey size ratio) and spatiotemporal co-occurrence between a predator and its prey (reflected in the accessibility coefficients, Fu et al., 2013). The minimum and maximum predator/prey size ratios were computed from the literature when available, or derived from observed diets and species' mean sizes, differing among size classes to account for ontogenetic changes in feeding behaviour (Supplementary Table S1). The amount of prey eaten by a given predator depends on the availability of prey. If enough prey items are present in a spatial cell, a predator feeds upon them uniformly until it reaches satiation, determined by the annual maximum ingestion rate taken from Ainsworth (2006). Otherwise, the predator depletes all prey but suffers from different degrees of starvation mortality depending on the prey availability (Shin and Cury, 2004). Fishing mortality is assumed to be knife-edged (i.e. vulnerable to fishing when schools reach the age of recruitment to the fishery) and spatially homogeneous within each management area. At the end of each time step, the reproduction process is modelled by adding eggs to the modelled system according to the spawning biomass, age of sexual maturity, sex ratio, the relative annual fecundity, and spawning seasonality (Fu et al., 2013).

Consequently, input information needed for the 10 populations includes: (i) growth, reproduction and mortality parameters (Fu et al., 2013; Table 2); (ii) feeding size ranges expressed as minimum and maximum predator/prey size ratios (Supplementary Table S1); (iii) annual fishing mortality time series obtained from stock assessments for Herring (DFO, 2014), Cod (Forrest et al., 2015), Lingcod (King et al., 2012) and from approximations based on catch data and biomass estimates for Arrowtooth Flounder and Walleye Pollock (Supplementary Table S2); (iv) fishing and reproduction seasonality (Supplementary Table S3); and (v) distribution maps for different life stages and time steps.

For the background taxa, only the predation, spatial distribution and movement processes are simulated. At the beginning of each year, the biomass of background taxa is separated into young-of-the-year, juveniles, and adults. The biomass in each life stage is then converted to abundance based on the average weight of each life stage and subsequently divided into 120 schools. At each time step, these schools at different life stages interact with schools of

<table><caption>**Table 2.** Growth, reproduction and mortality parameters for each of the key species considered in the PNCIMA OSMOSE model.</caption><thead><tr><th rowspan="2">Species</th><th colspan="5">Growth</th><th colspan="2">Reproduction</th><th colspan="3">Survivalship</th></tr><tr><th>L<sub>∞</sub> (cm)</th><th>K(year<sup>-1</sup>)</th><th>t<sub>0</sub> (year)</th><th>c(g·cm<sup>-3</sup>)</th><th>b</th><th>φ(eggs·g<sup>-1</sup>)</th><th>A<sub>mat</sub>(year)</th><th>A<sub>max</sub>(year)</th><th>A<sub>rec</sub>(year)</th><th>M<sup>diverse</sup> (year<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Euphausiids</td><td>1.84</td><td>1.68</td><td>-0.20</td><td>0.0091</td><td>2.920</td><td>24 469</td><td>0.3</td><td>1.7</td><td>0.5</td><td>0.485</td></tr><tr><td>Pacific Herring</td><td>26.3</td><td>0.36</td><td>-0.03</td><td>0.0070</td><td>2.997</td><td>200</td><td>3</td><td>10</td><td>3</td><td>0.130</td></tr><tr><td>Arrowtooth Flounder</td><td>58.92</td><td>0.28</td><td>0.48</td><td>0.0036</td><td>3.251</td><td>743</td><td>5</td><td>25</td><td>5</td><td>0.120</td></tr><tr><td>Walleye Pollock</td><td>44.50</td><td>0.92</td><td>0.57</td><td>0.0065</td><td>2.997</td><td>300</td><td>3</td><td>10</td><td>3</td><td>0.250</td></tr><tr><td>Pacific Cod</td><td>89.48</td><td>0.31</td><td>-0.12</td><td>0.0074</td><td>3.096</td><td>564</td><td>3</td><td>10</td><td>3</td><td>0.255</td></tr><tr><td>Lingcod</td><td>112.80</td><td>0.15</td><td>-3.01</td><td>0.0013</td><td>3.324</td><td>26</td><td>4.5</td><td>17</td><td>5</td><td>0.135</td></tr></tbody></table>

Growth parameters include $L_{∞}$, $k$, and $t_0$ for the von Bertalanffy growth model as well as parameters $c$ and $b$ for the weight-at-length allometric function. Relative fecundity $\phi$ is the number of eggs spawned per gram of mature female per year. $A_{mat}$, $A_{max}$, $A_{rec}$ are, respectively, the age at sexual maturity, the longevity, and the age of recruitment into the fisheries. The parameter $M^{diverse}$ is the mortality due to disease, senescence and predation by organisms not represented in the OSMOSE model (e.g. birds).
the 10 key populations and other background taxa through predator-prey relationships by exerting predation mortality or providing food. Mortality other than predation mortality, growth, and reproduction are not modelled for background taxa. Input information needed for the background taxa includes: (i) biomass time series from different sources, including stock assessments, survey estimates, and also estimates from Ainsworth (2006); (ii) mean length and weight at different stages; (iii) minimum and maximum predator/prey size ratios (Supplementary Table S1); and (iv) distribution maps for different life stages and time steps.

LTL groups only serve as food in PNCIMA OSMOSE and are represented as spatially distributed biomass pools. Input information for LTL groups includes: (i) biomass time series; (ii) minimum and maximum body sizes, and trophic levels (Supplementary Table S4); and (iii) distribution maps at different time steps derived from a number of different sources: http://www.science.oregonstate.edu/ocean productivit/; http://gdata1.sci.gsfc.nasa.gov/; Continuous Plankton Recorder data available at http://www.pices.int/projects/tcprsofp/data.aspx; and plankton databases (M. Galbraith, pers. comm.). LTL biomass observations covered short time periods and limited space; therefore, we used annual variability in the cumulative upwelling index (CUI) as a proxy for LTL biomass variations over time; Preikshot (2005) stated that the BC shelf primary production was more correlated with the upwelling index than with the Pacific Decadal Oscillation (PDO). CUI data were obtained from buoy observations (M. Foremen, pers. comm.). Based on CUI data, the Total Upwelling Magnitude Index (TUMI), a measure of the intensity of coastal upwelling integrated over the entire length of the defined upwelling season, was calculated according to Bograd et al. (2009) and then standardized to obtain time series of TUMI anomalies (TUMI_t). Annual LTL biomass time series B_t^LTL were then calculated as:

$$B_t^{LTL} = B_0^{LTL} \cdot e^{\delta \cdot \text{TUMI}_t},$$

where B_0^LTL is the 1950 LTL biomass level from Ainsworth (2006), and δ is a parameter that determines the degree of impact from TUMI.

## Simulations

The M^diverse of the first life stage (eggs and first-feeding larvae), referred to as M^diverse, is due to different causes (e.g. non-fertilization of eggs, starvation of first feeding larvae, advection, sinking, and predation by organisms not represented in OSMOSE) and is usually very hard to quantify. The first step of running the PNCIMA OSMOSE model was to estimate the M^diverse values for the 10 populations, as in Fu et al. (2013), so that the simulated biomass time series of the populations were as close as possible to those from stock assessments (Herring, Cod, and Lingcod), or when stock assessments were not available (Euphausiids, Arrowtooth Flounder, and Walleye Pollock), as close as possible to the biomass reported in Ainsworth (2006). This procedure is similar to the calibration of Ecosim models, where the vulnerability parameters of functional groups are estimated based on an iterative search (Monte-Carlo simulations) so that the simulated biomass and catch time series are as close as possible to the observed (Christensen and Walters, 2004). Although optimization methods have been developed for Osmose applications (Duboz et al., 2010; Oliveros-Ramos and Shin, 2016; Oliveros-Ramos

et al., 2017), these have not been applied here yet because of versions compatibility. In this article, we present results on Herring and Cod for which stock assessments have been conducted in the last 5 years.

Each PNCIMA OSMOSE simulation was run 100 times from 1940 to 2014 and outputs were averaged over the 100 realizations. The period 1940–1950 was the burn-in period. For the 10 populations, rather than biomass time series, only initial biomass was required to run the model. The initial biomasses of Arrowtooth Flounder, Walleye Pollock, and Euphausiids were the levels defined for year 1950 in Ainsworth (2006). Herring populations were initialized at the biomass levels of 1951, the first year for which stock assessments became available. Because Cod stock assessments began in 1956, PNCIMA OSMOSE initialized Cod populations at the virgin biomass level estimated from stock assessments. The initial biomass levels of Lingcod populations in year 1940 were taken directly from stock assessment.

In order to further improve the fittings to the biomass time series from stock assessments, we incorporated annual variability in LTL biomass as stated above. Furthermore, understanding that M^diverse is non-static and subject to variability in environmental conditions, such as upwelling variability (Ware and Thomson, 1991), we hypothesized that the M^diverse values of Cod and Herring populations were correlated with the TUMI anomalies:

$$M_{0,t}^{\text{diverse}} = M_0^{\text{diverse}} \cdot e^{\delta \cdot \text{TUMI}_t},$$

thus accommodating temporal variability in this parameter.

In order to find the optimal δ values to reflect the degree of variability in LTL biomass, and larval mortality of Cod and Herring, we ran the OSMOSE simulations using different δ values. Each δ value resulted in a series of residual (i.e. the difference between simulated and assessed biomass estimates) sum of squares (RSS) for each population of Herring and Cod. RSS values for each population under different δ values were scaled to obtain relative difference (RD):

$$RD = \frac{RSS_{\delta} - RSS_{\min}}{RSS_{\min}},$$

where RSS_δ is each individual RSS at a certain value of δ; RSS_min is the minimum value among all the RSS_δ values. The computation of RD allowed us to eliminate the effect of different scales in the RSS values among the different populations. The optimal δ value corresponds to the lowest sum of RD over all populations.

There were three steps to select an optimal δ value for LTL biomass, and larval mortality of Cod and Herring, respectively: (i) using a different δ value for LTL biomass, ranging from -0.5 to 0.5 (Table 3a) with a negative value implying an inverse relationship between LTL biomass and TUMI anomalies; (ii) keeping δ for LTL biomass at the optimal value obtained at Step 1 but varying δ for Cod larval mortality (Table 3b); (iii) keeping δ at the optimal values for LTL biomass and Cod larval mortality determined at Step 2 but varying δ for Herring larval mortality (Table 3c). After the three optimal δ values were determined, the resultant “Best” scenario was then used for further simulations.

Based on the Best scenario with the optimal δ values, we examined predation mortality time series for the five key populations (i.e. Herring-PRD, Herring-HG, Herring-CC, Cod-QCS, and Cod-HS) to identify temporal variability or trends. We then identified predators that were responsible for the predation mortality
<table><caption>**Table 3.** RSS and RD of RSS for each population of the three key species (Cod in QCS and HS; Herring in PRD, HG, and CC) and total RD summed over all populations (RD Sum) under three sets of simulations: a. using a different δ value for low-trophic-level (LTL) biomass variability; b. keeping δ for LTL biomass at 0.2 and varying δ for Cod larval mortality; c. keeping δ for LTL biomass at 0.2 and Cod larval mortality at 0.1 but varying δ for Herring larval mortality.</caption><thead><tr><th rowspan="3"></th><th colspan="10">RSS of Individual Population</th><th colspan="10">RD of Individual Population</th></tr><tr><th colspan="3">Cod</th><th colspan="3">Herring</th><th colspan="3">Cod</th><th colspan="3">Herring</th><th colspan="3">RD Sum</th></tr><tr><th>δ value</th><th>QCS</th><th>HS</th><th>PRD</th><th>HG</th><th>CC</th><th>QCS</th><th>HS</th><th>PRD</th><th>HG</th><th>CC</th><th>PRD</th><th>HG</th><th>CC</th><th>PRD</th><th>HG</th><th>CC</th><th>PRD</th><th>HG</th><th>CC</th><th>PRD</th><th>HG</th><th>CC</th></tr></thead><tbody><tr><td rowspan="15">a.</td><td>-0.5</td><td>10 044</td><td>50 104</td><td>194 997</td><td>150 528</td><td>347 226</td><td>12.85</td><td>7.58</td><td>1.11</td><td>0.63</td><td>4.51</td><td>26.68</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>-0.25</td><td>3103</td><td>20 156</td><td>134 979</td><td>98 504</td><td>161 046</td><td>3.28</td><td>2.45</td><td>0.46</td><td>0.07</td><td>1.56</td><td>7.82</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>-0.15</td><td>2494</td><td>14 252</td><td>120 696</td><td>93 738</td><td>114 397</td><td>2.44</td><td>1.44</td><td>0.31</td><td>0.01</td><td>0.82</td><td>5.02</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>-0.1</td><td>1457</td><td>9475</td><td>113 889</td><td>92 451</td><td>96 711</td><td>1.01</td><td>0.62</td><td>0.24</td><td>0</td><td>0.54</td><td>2.40</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>-0.05</td><td>2654</td><td>8617</td><td>104 601</td><td>94 794</td><td>83 985</td><td>2.66</td><td>0.48</td><td>0.13</td><td>0.03</td><td>0.33</td><td>3.63</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>0</td><td>784</td><td>8615</td><td>100 099</td><td>96 379</td><td>74 427</td><td>0.08</td><td>0.48</td><td>0.09</td><td>0.04</td><td>0.18</td><td>0.87</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>0.05</td><td>830</td><td>7519</td><td>99 134</td><td>99 153</td><td>68 422</td><td>0.14</td><td>0.29</td><td>0.08</td><td>0.07</td><td>0.09</td><td>0.67</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>0.1</td><td>1358</td><td>6418</td><td>92 198</td><td>102 412</td><td>64 479</td><td>0.87</td><td>0.10</td><td>0</td><td>0.11</td><td>0.02</td><td>1.10</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>0.15</td><td>867</td><td>6538</td><td>93 504</td><td>116 538</td><td>65 074</td><td>0.19</td><td>0.12</td><td>0.01</td><td>0.26</td><td>0.03</td><td>0.62</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td><strong>0.2</strong></td><td>725</td><td>6531</td><td>95 002</td><td>111 317</td><td>62 972</td><td>0.00</td><td>0.12</td><td>0.03</td><td>0.20</td><td>0</td><td><strong>0.35</strong></td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>0.25</td><td>893</td><td>5838</td><td>94 930</td><td>117 973</td><td>65 011</td><td>0.23</td><td>0</td><td>0.03</td><td>0.28</td><td>0.03</td><td>0.57</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>0.3</td><td>742</td><td>6295</td><td>100 447</td><td>121 175</td><td>68 301</td><td>0</td><td>0.08</td><td>0.09</td><td>0.31</td><td>0.08</td><td>0.59</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>0.35</td><td>795</td><td>6194</td><td>105 665</td><td>128 999</td><td>75 125</td><td>0.10</td><td>0.06</td><td>0.15</td><td>0.40</td><td>0.19</td><td>0.89</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>0.5</td><td>1040</td><td>6558</td><td>120 634</td><td>148 143</td><td>89 477</td><td>0.43</td><td>0.12</td><td>0.31</td><td>0.60</td><td>0.42</td><td>1.89</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td rowspan="10">b.</td><td>-0.2</td><td>1 915 682</td><td>909 961</td><td>66 629</td><td>116 841</td><td>72 216</td><td>2486.11</td><td>146.49</td><td>0.0002</td><td>0.03</td><td>0.12</td><td>2632.75</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>-0.15</td><td>714 636</td><td>397 211</td><td>66 618</td><td>116 505</td><td>71 721</td><td>926.81</td><td>63.38</td><td>0</td><td>0.02</td><td>0.11</td><td>990.32</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>-0.1</td><td>148 439</td><td>127 676</td><td>75 404</td><td>117 135</td><td>71 429</td><td>191.72</td><td>19.69</td><td>0.13</td><td>0.03</td><td>0.10</td><td>211.68</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>-0.05</td><td>10 718</td><td>25 270</td><td>85 312</td><td>117 460</td><td>70 555</td><td>12.91</td><td>3.10</td><td>0.28</td><td>0.03</td><td>0.09</td><td>16.41</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>-0.025</td><td>2104</td><td>13 217</td><td>88 672</td><td>117 957</td><td>66 844</td><td>1.73</td><td>1.14</td><td>0.33</td><td>0.03</td><td>0.03</td><td>3.27</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>0</td><td>1387</td><td>9477</td><td>96 717</td><td>116 952</td><td>66 706</td><td>0.80</td><td>0.54</td><td>0.45</td><td>0.03</td><td>0.03</td><td>1.85</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>0.025</td><td>816</td><td>8467</td><td>99 429</td><td>117 191</td><td>65 534</td><td>0.06</td><td>0.37</td><td>0.49</td><td>0.03</td><td>0.01</td><td>0.97</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>0.05</td><td>770</td><td>7582</td><td>98 154</td><td>116 812</td><td>67 624</td><td>0.00</td><td>0.23</td><td>0.47</td><td>0.02</td><td>0.05</td><td>0.77</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td><strong>0.1</strong></td><td>792</td><td>6170</td><td>93 803</td><td>116 163</td><td>66 648</td><td>0.03</td><td>0</td><td>0.41</td><td>0.02</td><td>0.03</td><td><strong>0.49</strong></td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>0.15</td><td>35 771</td><td>32 390</td><td>91 375</td><td>113 969</td><td>64 684</td><td>45.44</td><td>4.25</td><td>0.37</td><td>0</td><td>0</td><td>50.06</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>0.2</td><td>624 940</td><td>220 546</td><td>84 317</td><td>115 767</td><td>65 047</td><td>810.35</td><td>34.75</td><td>0.27</td><td>0.02</td><td>0.01</td><td>845.39</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td rowspan="10">c.</td><td>-0.15</td><td>879</td><td>6345</td><td>151 392</td><td>142 623</td><td>118 057</td><td>0.20</td><td>0</td><td>0.65</td><td>0.37</td><td>0.94</td><td>2.154</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>-0.1</td><td>731</td><td>6373</td><td>127 771</td><td>129 378</td><td>92 321</td><td>0.00</td><td>0.00</td><td>0.39</td><td>0.24</td><td>0.52</td><td>1.150</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>-0.05</td><td>936</td><td>6390</td><td>102 335</td><td>119 156</td><td>75 668</td><td>0.28</td><td>0.01</td><td>0.11</td><td>0.14</td><td>0.24</td><td>0.784</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>-0.025</td><td>830</td><td>6747</td><td>97 607</td><td>114 188</td><td>69 448</td><td>0.13</td><td>0.06</td><td>0.06</td><td>0.09</td><td>0.14</td><td>0.495</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>0</td><td>925</td><td>6531</td><td>95 002</td><td>111 317</td><td>62 972</td><td>0.27</td><td>0.03</td><td>0.03</td><td>0.07</td><td>0.03</td><td>0.429</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td><strong>0.025</strong></td><td>767</td><td>6407</td><td>91 983</td><td>107 178</td><td>60 849</td><td>0.05</td><td>0.01</td><td>0</td><td>0.03</td><td>0</td><td><strong>0.085</strong></td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>0.05</td><td>1361</td><td>6908</td><td>94 595</td><td>105 308</td><td>62 197</td><td>0.86</td><td>0.09</td><td>0.03</td><td>0.01</td><td>0.02</td><td>1.009</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>0.1</td><td>1053</td><td>7686</td><td>95 375</td><td>105 207</td><td>68 295</td><td>0.44</td><td>0.21</td><td>0.04</td><td>0.01</td><td>0.12</td><td>0.818</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>0.15</td><td>894</td><td>8046</td><td>104 024</td><td>104 386</td><td>75 466</td><td>0.22</td><td>0.27</td><td>0.13</td><td>0.00</td><td>0.24</td><td>0.862</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr></tbody></table>

The lowest RD Sum and its corresponding δ value are in bold.

of these five populations and assessed how the predation pressure from these predators changed over time. Predation pressure on prey by a particular predator was measured as the percentage of the prey biomass consumed by this predator over the total prey biomass consumed by all the species and taxa explicitly considered in PNCIMA OSMOSE.

## Results

### Simulated biomass trajectories

The RD for each population under each δ value (Table 3) indicates how far the simulated biomass trajectory is from the assessed one. The simulated biomass trajectory of Cod appeared to be most sensitive to the δ value for Cod larval mortality. A δ value > 0.15 or < 0 caused the RD of Cod populations to increase considerably (Table 3b), meaning their biomass trajectories departed drastically from the assessed biomass time series. The Best scenario was

chosen with δ values being at 0.2 for LTL biomass (Table 3a), 0.1 for Cod larval mortality (Table 3b), and 0.025 for Herring larval mortality (Table 3c). For clearer illustration, we displayed the simulated biomass trajectories, along with the assessed ones, under the Best scenario and the scenarios with no-environmental forcing (No-Env scenario), i.e. δ value being at 0 for both LTL biomass and larval mortality. For Herring, simulated biomasses from the Best scenario were not very different from those under the No-Env scenario (Figure 3a-c). In contrast, the simulated Cod biomasses from the Best scenario followed the cycles in the assessed biomasses much better than the No-Env scenario (Figure 3d and e).

### Predation mortality and predation pressure under the Best scenario

Temporal patterns of predation mortality varied among the different populations. Predation mortality for Herring did not have
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Biomass estimates from stock assessment (Assessed), simulated biomass trajectories averaged over 100 runs of the PNCIMA OSMOSE model under the Best scenario (Best) as well as the scenario where no temporal environmental variability was imposed on LTL biomass and larval mortality (No-Env) for the Herring populations in (a) PRD, (b) HG, and (c) CC, and for the Cod populations in (d) QCS, and (e) HS.</figcaption></figure>

a temporal trend in any of the three populations (Figure 4a–c). However, the predation mortality for Herring-PRD was consistently smaller than that for Herring-HG and Herring-CC. The first two decades of simulation were marked with predation mortality peaks for Herring-HG population. In contrast to Herring populations, the Cod populations displayed increasing trends in their predation mortality (Figure 4d and e). Moreover, the predation mortality for Cod-QCS was higher than that for Cod-HS.

Predation pressure from a specific predator species/taxon changed over time for all the populations of Herring and Cod. For clarity purposes, we averaged the predation pressure for three periods: (i) 1951–1966, the initial period when Herring populations were under heavy fishing mortality, marine mammals were

at low population levels, and Cod biomass remained high; (ii) 1983–1992, when Cod biomass experienced another rebound and marine mammals started to recover; and (iii) 2004–2014, the last decade of simulation, during which fishing mortality for Herring was minimal and predation mortality was relatively high for Cod.

Under the Best scenario, during the period 1951–1966, Herring-PRD suffered the highest predation pressure from Cod-HS (31%), which had subsequently decreased during the latter two periods (Figure 5a–c). Similarly, predation pressure on Herring-PRD from Lingcod-HS had also progressively reduced from 24 to 9%. In contrast, predation pressure on Herring-PRD from seals and sea lions progressively increased from 5% during the first period to 27% during the third period. Predation
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Trajectories of predation mortality (year⁻¹) averaged over 100 runs of the PNCIMA OSMOSE model for the Herring populations in (a) PRD, (b) HG, and (c) CC, and for the Cod populations in (d) QCS, and (e) HS.</figcaption></figure>

pressure from whales also increased drastically from 0% during the first period to 18% during the third period. Arrowtooth Flounder was an important predator consistently throughout the three periods with a narrow range from 9 to 12%.

For Herring-HG, Pacific Hake (hereafter abbreviated as Hake) (*Merluccius productus*) was one of the most important predators during all three periods of the simulation (Figure 5d–f). Predation pressure from marine mammals (whales, seals, and sea lions) on Herring-HG was not as dominant as on Herring-PRD, but its increasing trend over the three periods was consistent with that of Herring-PRD. Aside from marine mammals and Hake, Cod-HS, Lingcod-HS, and Spiny Dogfish (*Squalus acanthias*) were also important predators of Herring-HG in the Best scenario.

Hake became even more important as a predator of Herring moving southwards. Hake exerted 42% of the predation pressure on Herring-CC during the first period of the simulation although this value reduced to 24% during the third period (Figure 5g–i). Aside from Hake, marine mammals were the most important

predators of Herring-CC, accounting for 49% of the predation pressure during the third period of the simulation, increasing from 7% during the first period.

For Cod-QCS, Hake was the most important predator during the first period of simulation but its predation pressure decreased over the three periods (Figure 6a–c). In contrast, the predation pressure from seals and sea lions progressively increased to 46% during the third period. Similarly, the predation pressure from seals and sea lions on Cod-HS increased from 15% during the first period to 42% during the third period (Figure 6d–f). Lingcod was the most dominant predator for Cod-HS during the first period of simulation (30%, Figure 6d), however its dominance reduced in the latter two periods (Figure 6e and f). In addition, Pacific Halibut (hereafter abbreviated as Halibut) (*Hippoglossus stenolepis*) was also an important predator on Cod during the second period, but it became much less influential during the third period.

Overall, the PNCIMA OSMOSE model provided us with good pictures of (i) the levels of predation mortality for the different populations of the modelled key species and (ii) species and taxa that contribute to the predation mortality of different populations displayed spatial heterogeneity.

# Discussion

## Ecosystem model and validation

Through the development of the PNCIMA OSMOSE ecosystem model, we were able to reconstruct the population dynamics of all modelled species in different management areas, examine factors influencing their dynamics including species interactions in the form of predation mortality and environmental forcing on LTL biomass and larval mortality, and provide modelling and data/knowledge bases for Management Strategy Evaluation in an ecosystem context (Grüss et al., 2016b).

Biomass estimates from single-species-based stock assessment models represent the state-of-the-art in our understanding of the population status and have been used for setting harvest levels in fisheries management despite the potentially large estimation uncertainties. In the case of Herring and Cod, we used biomass estimates from the most recent stock assessments for validating the PNCIMA OSMOSE ecosystem model. If more accurate data or estimates become available, we will further validate our ecosystem model against the new information. Through the validation process, we identified optimal ways to incorporate temporal variability into LTL biomass and larval mortality via the TUMI anomalies, and in particular, we were able to reproduce the cycles in the assessed biomasses (Forrest et al. 2015). In the same manner (i.e. validation testing), other environmental factors (e.g. PDO) and forms of functionality can be investigated to promote a more mechanistic understanding of the ecosystem and to further improve the PNCIMA OSMOSE model.

The PNCIMA OSMOSE model generated good fitting to the assessed biomass time series for Herring and Cod, nevertheless, there were notable discrepancies. Aside from the possibility that single-species-based assessments do not account for some processes that an ecosystem assessment would, these discrepancies could be due to a number of reasons. Firstly, for Herring, preliminary analyses of historical length-at-age data suggested that the von Bertalanffy growth parameters had changed over time (DFO, 2014), likely caused by changes in temperature. If such temporal changes are imposed on the growth parameters of the PNCIMA
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Contributors to the predation mortality of Herring predicted by the PNCIMA OSMOSE model for the PRD area (a–c), the HG area (d–f) and the CC area (g, h, i). Panels (a, d, and g) are for the period from 1951 to 1966; panels (b, e, and h) for the period from 1983 to 1992; and panels (c, f, and i) for the period from 2004 to 2014. Key species/background taxa contributing to < 5% of the predation mortality of Herring in all of the three Herring management areas and during all three periods were combined into the Other Group (Oth) with species including Walleye Pollock, Petrale Sole, Slope Rockfish, Inshore Rockfish, Sablefish, and Coho and Chinook Salmon. One hundred simulation runs were used to produce pie charts. SS, Seals and sea lions; Wha, Whales; L-HS, Lingcod in HS; C-HS, Cod-HS; L-QCS, Lingcod in QCS; AF, Arrowtooth Flounder; ShR, Shelf Rockfish; Hal, Halibut; DF, Spiny Dogfish.</figcaption></figure>

OSMOSE model, the fitting may be improved. Alternatively, the effect of temperature change can be modelled at the individual physiological level (e.g. Vasseur and McCann, 2005). Secondly, although the fitting for Herring has been improved by correlating LTL biomasses with the TUMI anomalies, information on the initiation of upwelling, phytoplankton bloom, and food availability to larvae due to spatiotemporal match and mismatch is still unavailable for further improvement. It is worth noting that in the Strait of Georgia, adjacent to PNCIMA, the timing of the copepod biomass peak has shifted unequivocally earlier every year, which has caused reproductive failure among some marine bird species and may have affected the survival of juvenile fish (Johannessen and McCarter, 2010). In PNCIMA, climate change has produced shifts in dominant copepod species that altered food nutritional quality for juvenile fish (Cummins and Haigh, 2010). Such climate-driven changes in the plankton community should be monitored continually; their potential impacts on

higher trophic levels, currently not known (Johannessen and McCarter, 2010), can be investigated through ecosystem models such as OSMOSE. Thirdly, the assumption of homogeneous fishing mortality over space might limit model performance, particularly in the case of Herring, because fishing fleets tend to target Herring spawning aggregations, which might have effectively removed a large proportion of the spawning biomass and contributed to the collapse of Herring stocks (Council of the Haida Nation, 2011). Model development to introduce spatial heterogeneity of fishing mortality in OSMOSE is underway, which will allow the dynamics of fishing fleets to be adequately simulated through the use of distribution maps of fishing fleets and effort.

## Spatial population structure and background taxa

Through modelling the different spatially structured populations, we were able to illustrate differences in predation mortality and
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Contributors to the predation mortality of Cod predicted by the PNCIMA OSMOSE model for the QCS area (a–c) and the HS area (d–f). Panels (a and d) are for the period from 1951 to 1966; panels (b and e) for the period from 1983 to 1992; and panels (c and f) for the period from 2004 to 2014. Key species/background taxa contributing to < 5% of the predation mortality of Cod in all of the two management areas and during all three period were combined into the Oth with species including Shelf Rockfish, Arrowtooth Flounder, Walleye Pollock, Petrale Sole, Slope Rockfish, Inshore Rockfish, and Sablefish. One hundred simulation runs were used to produce pie charts. SS, Seals and sea lions; Wha, Whales; Hal, Halibut; DF, Spiny Dogfish. Ling and Cod in panels (a–c) represent Lingcod-QCS and Cod-QCS; Ling and Cod in panels (d–f) represent Lingcod-HS and Cod-HS.</figcaption></figure>

contributing predators among the different populations of the same species. For example, predation mortality for Herring-PRD was generally lower than that for Herring-HG and Herring-CC over the period 1951–2014. If spatial heterogeneity in predation mortality was not acknowledged and the populations were managed as one population, fishing levels for Herring in different spatial regions might be inappropriately set. Similarly, for Cod-QCS, the average predation mortality was significantly higher and the increasing trend in predation mortality was also more pronounced than those for Cod-HS. This result, combined with the fact that historical tagging data indicated that Cod rarely moved between QCS and HS (Forrest et al., 2015), justifies the assessment and management of Cod-QCS and Cod-HS on an area-specific basis.

Including 19 background taxa in the PNCIMA OSMOSE model allowed us to explicitly account for all the known potential predators of the key species, thus permitting a comprehensive evaluation of the predation pressure on the key species in PNCIMA. For example, the PNCIMA OSMOSE model allowed us to acknowledge that non-resident taxa such as Hake and Humpback Whale (Megaptera novaeangliae), whose entire life cycle would be difficult to model with OSMOSE, were responsible for a high proportion of the predation mortality of Herring populations, in agreement with Ford et al. (2009) and Schweigert et al. (2010) who drew empirical conclusions based on direct

abundance observations of studied species. In future versions of the PNCIMA OSMOSE model we will consider the fishing mortality and average growth of the background taxa to provide a more accurate picture of ecosystem dynamics.

The contributions of different species/taxa to the predation mortality of Herring estimated in this study and those reported in Schweigert et al. (2010) were generally in concordance, although the two studies considered two different but adjacent geographical areas. For instance, both studies indicated that marine mammals had become increasingly important predators of Herring in the most recent years. Schweigert et al. (2010) found that Halibut consumed the least amount of Herring off the West Coast of Vancouver Island (WCVI). In contrast, the PNCIMA OSMOSE simulations indicated that Halibut could be an important predator of Herring, although its importance varied over different areas resulting from different degrees of spatial overlap between these two species. The predation pressure exerted by Halibut on Herring over the entire period was lowest in CC (adjacent to WCVI), and was highest in HG, particularly during the period 1951–1966. Traditional local knowledge (Council of the Haida Nation, 2011) reported frequent observations of Halibut following Herring in HG for the purpose of feeding.

Schweigert et al. (2010) reported that 13.3% of Herring biomass was consumed by Lingcod in WCVI, which compared well with PNCIMA OSMOSE prediction for CC during the second
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Starvation mortality (year⁻¹) time series averaged over 100 runs of the PNCIMA OSMOSE model for (a) the Herring populations of PRD, HG, and CC; and (b) the Cod populations of HS and QCS.</figcaption></figure>

period (1983–1992; 13% with 10% by Lingcod-QCS and 3% by Lingcod-HS). In the northern regions of PNCIMA (i.e. PRD and HG), the predation pressure of Lingcod-HS on Herring was larger, particularly during the first two periods. Similarly, Cod was a more important predator of Herring in the northern regions of PNCIMA than in CC.

Schweigert et al. (2010) found that 9.6% of the biomass of Herring was consumed by Hake in WCVI, a figure that was much smaller than OSMOSE predicted for CC, which could be due to the fact that the biomass estimates of Hake were based on biomass indices rather than total biomass, representing minimum estimates (Schweigert et al., 2010). In contrast, the biomass time series of Hake used in the present study were based on recent stock assessments (International Joint Technical Committee for Pacific Hake, 2013).

Another study examining stomach contents from HS indicated that 49% of the diet (by biomass) of Petrale Sole (*Eopsetta jordani*) consisted of Herring (Pearsall and Fargo, 2007). In PNCIMA OSMOSE, the consumption of Herring by Petrale Sole in all management areas was minimal due to the small biomass of Petrale Sole. In other words, although Herring is important as a food source, at its current population level Petrale Sole does not impose a threat on Herring.

## Factors influencing population dynamics

The aforementioned knowledge of predation mortality and contributing predators obtained from the PNCIMA OSMOSE model can provide us with insights on the dynamics of the key populations. Herring populations in PNCIMA have been depressed during the past decade despite annual fishing closures in HG and a drastic reduction in fishing activities in PRD and CC. It was suggested that ecosystem forcing, such as food supply (bottom-up), predation (top-down), or competition may affect Herring populations more than fishing (Schweigert et al., 2010). Herring stock assessments indicated that the natural mortality of Herring had been increasing in PRD, HG, and CC in the last two decades

(DFO, 2014). The PNCIMA OSMOSE simulations indicated that predation mortality of the different Herring populations did not have increasing trends; however, starvation mortality displayed increasing trends (Figure 7a), which suggested that food availability could be a limiting factor for the Herring populations.

Cod stock assessment, assuming constant natural mortality over time but estimating time-varying age-2 recruits and annual fishing mortality, indicated that age-2 recruits declined over time, but was unable to identify the reasons (Forrest et al., 2015). The PNCIMA OSMOSE simulations suggested that predation mortality, primarily on pre-recruits, increased over time significantly for both Cod-QCS and Cod-HS populations. Starvation mortality on the other hand was negligible and did not show an increasing trend (Figure 7b), suggesting increased predation mortality rather than food availability may have been limiting the Cod populations.

The PNCIMA OSMOSE model demonstrated that the increases of marine mammals including whales, particularly seals and sea lions were responsible for the increased predation mortality for Cod but didn't afflict increased predation mortality for Herring, inconsistent with the hypothesis that the increase of marine mammals impeded the recovery of Herring populations (e.g. Schweigert et al., 2010; Surma and Pitcher, 2015). This was because the increase of marine mammals resulted in the reduction of other important predators of Herring, such as Cod, Lingcod, and Halibut through predation, which subsequently stabilized the predation mortality of Herring. Overall, modelling the PNCIMA ecosystem in its entirety by explicitly considering nearly all the species/taxa has presented to us a more complete picture of species interactions and their impacts on the ecosystem.

With the new OSMOSE modelling platform we will explore various scenarios of fishing strategies and environmental changes, similar to what was done in Fu et al. (2013). As an initial step, we chose four hypothetical scenarios (no fishing for Herring, fishing but no predation for Herring, no fishing for Cod, and fishing but no predation for Cod) to investigate the relative impacts of fishing and predation mortality on the dynamics of Herring and Cod
populations (Supplementary Appendix S1). Results indicated that the relative impacts of different factors (e.g. fishing and preda- tion) may depend on species, as well as populations of the same species. Future simulations involving more realistic scenarios will allow us to identify ecosystem-level factors impacting multiple species and fisheries over the long term that could be missed un- der single-species-based research.

## Summary

In conclusion, the enhanced version of the OSMOSE model has allowed us to explicitly investigate the role of predation on the dynamics of key populations of the PNCIMA ecosystem and to evaluate the contribution of the different species of this ecosystem to the predation mortality of the key populations over different time periods and for different spatial regions. These investigations provide invaluable insights into species interactions and popula- tion dynamics in the PNCIMA ecosystem that will help improve fisheries management off western Canada in an ecosystem context.

## Supplementary data

Supplementary material is available at the *ICESJMS* online version of the article.

## Acknowledgements

This research is funded by the Strategic Program for Ecosystem- Based Research and Advice in Fisheries and Oceans Canada. We are extremely grateful to the following colleagues (listed alphabet- ically) who have generously provided different types of data: Leslie Barton, Jennifer Boldt, Jaclyn Cleary, Kristen Daniel, John Ford, Robyn Forrest, Moira Galbraith, Chris Grandin, Kendra Holt, Jim Irvine, Sheena Majewski, John Morris, R. Ian Perry, Ian Stewart, Matt Thompson, and Marc Trudel. We are very grateful to Jason S. Link (Editor) and two anonymous journal peer reviewers, as well as Lisa Christensen (Science Advisor at Canadian Science Advisory Secretariat) for providing helpful comments and edits. Y.J.S. and P.V. were supported by the French project EMIBIOS (End-to-end Modelling and Indicators for BIOdiversity Scenarios) with Foundation pour la Recherche sur la Biodiversité (FRB, contract no. APP-SCEN-2010-II).

## References

* Ainsworth, C. 2006. Strategic marine ecosystem restoration in north- ern British Columbia. PhD thesis. Department of Resource Management and Environmental Studies. University of British Columbia, Vancouver, BC 422 pp.

* Beacham, T. D., Schweigert, J. F., MacConnachie, C., Le, K. D., and Flostrand, L. 2008. Use of microsatellites to determine population structure and migration of Pacific Herring in British Columbia and Adjacent Regions. Transactions of the American Fisheries Society, 137: 1795-1811.

* Bograd, S. J., Schroeder, I., Sarkar, N., Qiu, X., Sydeman, W. J., and Schwing, F. B. 2009. Phenology of coastal upwelling in the California Current. Geophysical Research Letters, 36: L01602.

* Christensen, V., and Walters, C. J. 2004. Ecopath with Ecosim: meth- ods, capabilities and limitations. Ecological Modelling, 172: 109-139.

* Council of the Haida Nation. 2011. Haida Marine Traditional Knowledge Study. Vol. 1-3. 90 + 107 + 122 pp.

* Cummins, P., and Haigh, R. 2010. Ecosystem status and trends report for north coast and hecate strait ecozone. DFO Canadian Science Advisory Secretariat Research Document, 2010/046. vi + 61 p.

* DFO. 2014. Stock assessment and management advice for british columbia pacific herring: 2014 status and 2015 forecast. DFO Canadian Science Advisory Secretariat Science Advisory Report, 2014/060. 21 pp.

* Duboz, R., Versmisse, D., Travers, M., Ramat, E., and Shin, Y. J. 2010. Application of an evolutionary algorithm to the inverse pa- rameter estimation of an individual-based model. Ecological Modelling, 221: 840-849.

* Flostrand, L. A., Schweigert, J. F., Daniel, K. S., and Cleary, J. S. 2009. Measuring and modelling Pacific herring spawning-site fidelity and dispersal using tag-recovery dispersal curves. ICES Journal of Marine Science, 66: 1754-1761.

* Ford, J. K. B., Rambeau, A. L., Abernethy, R. M., Boogaards, M. D., Nichol, L. M., and Spaven, L. D. 2009. An Assessment of the Potential for Recovery of Humpback Whales off the Pacific Coast of Canada. DFO Canadian Science Advisory Secretariat Research Document, 2009/015. iv + 33 pp.

* Forrest, R. E., Rutherford, K. L., Lacko, L., Kronlund, A. R., Starr, P. J., and McClelland, E. K. 2015. Assessment of Pacific Cod ( Gadus macrocephalus ) for Hecate Strait (5CD) and Queen Charlotte Sound (5AB) in 2013. DFO Canadian Science Advisory Secretariat Research Document, 197. 2015/052. xii + pp.

* Fu, C., Perry, R. I., Shin, Y.-J., Schweigert, J., and Liu, H. 2013. An ecosystem modelling framework for incorporating climate regime shifts into fisheries management. Progress in Oceanography, 115: 53-64.

* Grimm, V., Ayllón, D., and Railsback, S. F. 2017. Next-generation in- dividual-based models integrate biodiversity and ecosystems: yes we can, and yes we must. Ecosystems, 20: 229-236.

* Grüss, A., Schirripa, M. J., Chagaris, D., Velez, L., Shin, Y.-J., Verley, P., Oliveros-Ramos, R., and Ainsworth, C. H. 2016a. Estimating natural mortality rates and simulating fishing scenarios for Gulf of Mexico red grouper ( Epinephelus morio ) using the ecosystem model OSMOSE-WFS. Journal of Marine Systems, 154: 264-279.

* Grüss, A., Harford, W. J., Schirripa, M. J., Velez, L., Sagarese, S. R., Shin, Y.-J., and Verley, P. 2016b. Management strategy evaluation using the individual-based, multi-species modeling approach OSMOSE. Ecological Modelling, 340: 86-105.

* International Joint Technical Committee for Pacific Hake. 2013. Status of the Pacific Hake (whiting) stock in U.S. and Canadian waters in 2013. 190 pp.

* Johannessen, S. C., and McCarter, B. 2010. Ecosystem status and trends report for the Strait of Georgia Ecozone. DFO Canadian Science Advisory Secretariat Research Document, 2010/010. vi + 43 pp.

* King, J. R., McAllister, M., Holt, K. R., and Starr, P. J. 2012. Lingcod ( Ophiodon elongatus ) stock assessement and yield advice for out- side stocks in British Columbia. DFO Canadian Science Advisory Secretariat Research Document, 2011/124. viii + 177 p.

* Lucas, B. G., Verrin, S., and Brown, R. (Eds) 2007. Ecosystem over- view: Pacific North Coast Integrated Management Area (PNCIMA). Canadian Technical Report of Fisheries and Aquatic Sciences, 2667. xiii + 104 pp.

* Oliveros-Ramos, R., and Shin, Y.-J. 2016. calibraR: an R package for the calibration of individual based models. arXiv:1603.03141.

* Oliveros-Ramos, R., Verley, P., Echevin, V., and Shin, Y. J. 2017. A sequential approach to calibrate ecosystem models with multiple time series data. Progress in Oceanography, 151: 227-244.

* Pearsall, I. A., and Fargo, J. J. 2007. Diet composition and habitat fi- delity for groundfish assemblages in Hecate Strait, British Columbia. Canadian Technical Report of Fisheries and Aquatic Sciences, 2692. vi + 141 pp.

* Preikshot, D. 2005. Data sources and derivation of parameters for generalised Northeast Pacific Ocean Ecopath with Ecosim models. University of British Columbia Fisheries Centre Research Reports, 13: 179-237 pp.
Schweigert, J. F., Boldt, J. L., Flostrand, L., and Cleary, J. S. 2010. A review of factors limiting recovery of Pacific Herring stocks in Canada. ICES Journal of Marine Science, 67: 1903–1913.

Seidl, R. 2017. To model or not to model, that is no longer the question for ecologists. Ecosystems, 20: 222–228.

Shin, Y.-J., and Cury, P. 2004. Using an individual-based model of fish assemblages to study the response of size spectra to changes in fishing. Canadian Journal of Fisheries and Aquatic Sciences, 61: 414–431.

Surma, S., and Pitcher, T. J. 2015. Predicting the effects of whale population recovery on Northeast Pacific food webs and fisheries: an ecosystem modelling approach. Fisheries Oceanography, 24: 291–305.

Travers-Trolet, M., Shin, Y.-J., and Field, J. G. 2014. An end-to-end coupled model ROMS-N2P2Z2D2-OSMOSE of the southern

Benguela foodweb: parameterisation, calibration and pattern-oriented validation. African Journal of Marine Science, 36: 11–29.

Travis, J., Coleman, F. C., Auster, P. J., Cury, P. M., Estes, J. A., Orensanz, J., Peterson, C. H., et al. 2014. Integrating the invisible fabric of nature into fisheries management. Proceedings of the National Academy of Sciences of the United States of America, 111: 581–584.

Vasseur, D. A., and McCann, K. S. 2005. A mechanistic approach for modeling temperature-dependent consumer-resource dynamics. The American Naturalist, 166: 184–198.

Ware, D. M., and Thomson, R. E. 1991. Link between long-term variability in upwelling and fish production in the northeast Pacific Ocean. Canadian Journal of Fisheries and Aquatic Sciences, 48: 2296–2306.

Handling editor: Jason Link