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

Abstract

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.

ecosystem modelling, fisheries management, OSMOSE, predation mortality, spatial population structure, species interaction

Introduction

Quantitative ecosystems models, i.e. mathematical representations of ecosystem dynamics, are valuable tools for studying complex interactions and their effects (Seidl, 2017). Species interactions and their spatial patterns are fundamental to the dynamics of fish populations and ecosystems, though often ignored in 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).

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).

We focused on six key species in the PNCIMA OSMOSE model: Pacific Herring (Clupea pallasii), Pacific Cod (Gadus 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.

Table 1.

Key species, background taxa and LTL groups included in the PNCIMA OSMOSE model.

Species/Taxon name	Type	Species represented
Euphausiids	Key species	Thysanoessa spp., Euphausia spp.
Pacific Herring	Key species	Pacific Herring (Clupea pallasi)
Arrowtooth Flounder	Key species	Arrowtooth Flounder (A. stomias)
Walleye Pollock	Key species	Walleye Pollock (T. chalcogramma)
Pacific Cod	Key species	Pacific Cod (G. macrocephalus)
Lingcod	Key species	Lingcod (O. elongatus)
Seals and sea lions	Background taxa	Harbour Seal (Phoca vitulina), Steller Sea Lion (Eumetopias jubatus)
Whales	Background taxa	Humpback Whale (M. novaeangliae)
Pacific Hake	Background taxa	Pacific Hake (M. productus)
Pacific Ocean Perch	Background taxa	Pacific Ocean Perch (Sebastes alutus)
Spiny Dogfish	Background taxa	Spiny Dogfish (S. acanthias)
Flatfish	Background taxa	Dover Sole (Microstomus pacificus), Rock Sole (Lepidosetta bilineata), English Sole (Parophyrys 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)
Pacific Halibut	Background taxa	Pacific Halibut (H. stenolepis)
Petrale Sole	Background taxa	Petrale Sole (E. jordani)
Shelf Rockfish	Background taxa	Yellowtail (Sebastes flavidus), Silvergray (Sebastes brevispinis), Widow (Sebastes entomelas), Bocaccio (Sebastes paucispinis), Canary (Sebastes pinniger)
Slope Rockfish	Background taxa	Yellowmouth (Sebastes reedi), Rougheye (Sebastes aleutioanus), Redstripe (Sebastes proriger), Sharpchin (Sebastes zacentrus), Redbanded (Sebastes babcocki), Shortspine Thornyhead (Sebastolobus altivelis), Splitnose (Sebastes diploproa), Longspine Thornyhead (Sebastolubus alascanus), Darkblotched (Sebastes cremeri), and Shortraker (Sebastes borealis)
Inshore Rockfish	Background taxa	Yelloweye (Sebastes rubberrimus), Quillback (Sebastes maliger), Copper (Sebastes caurinus), China (Sebastes nebulosus), Tiger (Sebastes nigrocinctus)
Spotted Ratfish	Background taxa	Spotted Ratfish (Hydrolagus colliei)
Sablefish	Background taxa	Sablefish (Anoplopoma fimbria)
Coho Chinook	Background taxa	Coho Salmon (Oncorhynchus kisutch), Chinook Salmon (Oncorhynchus tshawytscha)
Shallow Benthic Fish	Background taxa	Eelpouts (Zoarcidae), Poachers (Agonidae), Sculpins (Cottidae)
Forage Fish	Background taxa	Eulachon (Thaleichthys pacificus), Smelts (Osmeridae), Sandlance (Ammodytes hexapterus)
Crabs	Background taxa	Dungeness crab (Metacarcinus magister), Tanner Crab (Chionoecetes spp.), Red Rock Crab (Cancer productus)
Shrimp	Background taxa	Smooth Shrimp (Pandalus jordani), Spiny Shrimp (Pandalus borealis eous), Pink Shrimp (Pandalus goniturus), Sidestrip Shrimp (Pandalopsis disbar), Prawn (Pandalus platycterus)
Detritus Benthos	Background taxa
Phytoplankton	LTL group
Copepods	LTL group

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.

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.

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^{\mathrm{d}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{e}}$⁠) 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 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.

Table 2.

Growth, reproduction and mortality parameters for each of the key species considered in the PNCIMA OSMOSE model.

	Growth					Reproduction		Survivalship
Species	L∞ (cm)	K(year−1)	t0 (year)	c(g ⋅ cm−3)	b	φ(eggs ⋅ g−1)	Amat(year)	Amax(year)	Arec(year)	Mdiverse (year−1)
Euphausiids	1.84	1.68	−0.20	0.0091	2.920	24 469	0.3	1.7	0.5	0.485
Pacific Herring	26.3	0.36	−0.03	0.0070	2.997	200	3	10	3	0.130
Arrowtooth Flounder	58.92	0.28	0.48	0.0036	3.251	743	5	25	5	0.120
Walleye Pollock	44.50	0.92	0.57	0.0065	2.997	300	3	10	3	0.250
Pacific Cod	89.48	0.31	−0.12	0.0074	3.096	564	3	10	3	0.255
Lingcod	112.80	0.15	−3.01	0.0013	3.324	26	4.5	17	5	0.135

Growth parameters include $L_{\infty }$⁠, k, and t₀ for the von Bertalanffy growth model as well as parameters c and b for the weight-at-length allometric function. Relative fecundity $\varphi$ 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).

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 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.productivity/; http://gdata1.sci.gsfc.nasa.gov/; Continuous Plankton Recorder data available at http://www.pices.int/projects/tcprsotnp/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 (⁠${\mathrm{T}\mathrm{U}\mathrm{M}\mathrm{I}}_t$⁠). Annual LTL biomass time series $B_t^{\mathrm{L}\mathrm{T}\mathrm{L}}$ were then calculated as:

$$B_t^{\mathrm{L}\mathrm{T}\mathrm{L}}=\mathrm{ }{B}_0^{\mathrm{L}\mathrm{T}\mathrm{L}}·e^{\delta ·{\mathrm{T}\mathrm{U}\mathrm{M}\mathrm{I}}_t}$$

where $B_0^{\mathrm{L}\mathrm{T}\mathrm{L}}$ is the 1950 LTL biomass level from Ainsworth (2006), and $\delta$ is a parameter that determines the degree of impact from TUMI.

Simulations

The $M^{\mathrm{d}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{e}}$ of the first life stage (eggs and first-feeding larvae), referred to as $M_0^{\mathrm{d}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{e}}$⁠, 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_0^{\mathrm{d}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{e}}$ 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_0^{\mathrm{d}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{e}}$is non-static and subject to variability in environmental conditions, such as upwelling variability (Ware and Thornson, 1991), we hypothesized that the $M_0^{\mathrm{d}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{e}}$ values of Cod and Herring populations were correlated with the TUMI anomalies:

$${\mathrm{ }M}_{0,\mathrm{ }t}^{\mathrm{d}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{e}}=\mathrm{ }{M}_0^{\mathrm{d}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{e}}·e^{\delta ·{\mathrm{T}\mathrm{U}\mathrm{M}\mathrm{I}}_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 $\delta$ values. Each $\delta$ 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 $\delta$ values were scaled to obtain relative difference (RD):

$$RD=\mathrm{ }\frac{{RSS}_{\delta }-{RSS}_{\mathrm{m}\mathrm{i}\mathrm{n}}}{{RSS}_{\mathrm{m}\mathrm{i}\mathrm{n}}},$$

where ${RSS}_{\delta }$ is each individual RSS at a certain value of $\delta$⁠; ${RSS}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ is the minimum value among all the ${RSS}_{\delta }$ values. The computation of RD allowed us to eliminate the effect of different scales in the RSS values among the different populations. The optimal $\delta$ value corresponds to the lowest sum of RD over all populations.

There were three steps to select an optimal $\delta$ 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.

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.

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

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

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 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 $\delta$ 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 $\delta$ value for Cod larval mortality. A $\delta$ 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).

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.

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 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.

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.

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 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%.

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.

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.

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.

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 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 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 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.

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.

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 predation) 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 under 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 population 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 alphabetically) 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).