Ecosystem-based reference points under varying plankton productivity states and fisheries management strategies

Ecosystem-based reference points under varying plankton productivity states and ﬁsheries management strategies. In the context of ecosystem-based ﬁsheries management, which should consider changing and uncertain environmental conditions, the development of ecosystem-based biological reference points (EBRPs) to account for important multi-species (MS) interactions, ﬁshery operations, and climate change, is of paramount importance for sustainable ﬁsheries management. However, EBRPs under varying plankton productivity states and ﬁsheries management strategies are seldom developed, and the ecosystem effects of these changes are still largely unknown. In this study, ecosystem-based F MSY (ﬁshing mortality rate at MSY) values were estimated within an end-to-end ecosystem model (OSMOSE) for three focused ﬁsh species (Paciﬁc Herring, Clupea pallasii ; Paciﬁc Cod, Gadus macrocephalus ; Lingcod, Ophiodon elongatus ) under three plankton productivity states of differing plankton biomass at high, current, and low levels. In addition, ecosystem effects were compared across different plankton productivity and ﬁsheries management strategies with the latter consisting of two ﬁshery scenarios (i.e. single-species-focused (SS) and MS-focused), various ﬁshing mortality rates, and two harvest policies (with and without harvest control rules, HCRs). Main ﬁndings of this study include: (i) plankton productivity change affected the values of ecosystem-based F MSY , which increased as plankton productivity states changed from low to high plankton biomass; (ii) ecosystem-based F MSY for Paciﬁc Herring and Paciﬁc Cod stocks increased when ﬁshery scenarios shifted from SS-focused to MS-focused; (iii) ﬁsheries management incorporating HCR yielded more stable system catch and system biomass; and (iv) high plankton biomass combined with ﬁsheries management using HCR could maintain stable ecosystem production and sustainable ﬁsheries. Based on our ﬁndings, we highlight possible adaptive ﬁsheries management strategies in the face of future climate and ocean changes. Overall, EBRPs complement SS stock assessments by incorporating key ecological processes and ecosystem properties, thus providing supporting evidence for better incorporation of ecosystem considerations into scientiﬁc advice for sustainable ﬁsheries management.


Introduction
Biological reference points (BRPs) are a set of assigned or estimated objectives for a fish stock that represent biomass targets to achieve or thresholds for avoidance (Sainsbury, 2008;Punt et al., 2016). BRPs are generally used as benchmarks for indicating stock status and for identifying desirable levels of fishing mortality for achieving objectives. The most commonly used BRPs are based on the concept of MSY, i.e. the biomass or fishing mortality rate at MSY (i.e. B MSY , F MSY ) (Gabriel and Mace, 1999;Mace, 2001;Punt and Smith, 2001;Tyrrell et al., 2011). Over the last several decades, BRPs have been employed worldwide to support decisions for fisheries management (Braccini et al., 2015;Holsman et al., 2016). However, routinely used BRPs are mostly derived from SS equilibrium models, and do not account for ecosystemlevel ecological processes (e.g. species interactions, ecosystem productivity, or other environmental changes such as those arising from climate change). Failure to account for such uncertainties may reduce the effectiveness of applied BRPs (Gislason, 1999;Collie and Gislason, 2001;Holsman et al., 2016). Consequently, fisheries management decisions based on BRPs may become either more or less precautionary than originally intended (Heino et al., 2013). Therefore, the development and testing of ecosystem-based biological reference points (EBRPs) to account for important MS interactions (Smith et al., 2015;Dolan et al., 2016;Holsman et al., 2016), fishery operations (Forrest and Walters 2009;Forrest et al., 2015b), and climatic changes (Tyrrell et al., 2011;Moffitt et al., 2016;Payne et al., 2017) at the ecosystem level, is of paramount importance for moving towards sustainable ecosystem-based fisheries management (EBFM; Hall and Mainprize, 2004;Link, 2010) in the context of uncertain future changes. Development and testing of EBRPs is an important step in the move towards EBFM.
Despite increasing importance and awareness of the need for EBRPs, details on how EBRPs should be developed and assessed, how they differ under future scenarios of climate change and varying human activities (e.g. shifts in fishing strategies and harvest policies), and how ecosystems will be impacted by their applications are still not well-studied (Hall and Mainprize, 2004;Brunel et al., 2010;Link, 2010;Smith et al., 2015;Moffitt et al., 2016). However, there is increasing evidence that climate and ocean changes could have considerable impacts on marine ecosystems and fisheries sustainability (Brunel et al., 2010;Little et al., 2011;Sumaila et al., 2011;Heino et al., 2013;Smith et al., 2015).
Development and simulation-testing of EBRPs requires comprehensive and robust ecosystem modelling platforms (Plagányi, 2007;Fulton, 2010;Travers-Trolet et al., 2014). Ecosystem models can serve as excellent platforms for generating and evaluating EBRPs and quantifying the trade-offs between different management strategies in the face of climate change (Brown et al., 2010). Examples of such ecosystem models include Ecopath with Ecosim (Christensen and Walters, 2004), Atlantis (Fulton et al., 2005), and OSMOSE (Shin and Cury, 2004). Specifically, OSMOSE is an individual-based, spatially and temporally explicit modelling platform, able to account for size-based predator-prey interactions, whole life cycle dynamics of marine organisms, impacts of environmental changes on growth and survival of marine organisms, as well as the effects of harvest strategies and marine protected areas (Shin and Cury, 2004;Yemane et al., 2008;Travers-Trolet et al., 2014;Fu et al., 2017;Payne et al., 2017).
In the Northeast Pacific, climate and ocean changes has produced shifts in physical and chemical properties of the ocean, such as increased sea surface temperature, dissolved CO 2 and decreased dissolved oxygen Morrison et al., 2014), as well as altered lower-trophic-level (LTL) communities, resulting in alterations of food nutritional quality for juvenile fish (reviewed in Cummins et al., 2010). Plankton productivity, on which most marine animals depend directly or indirectly, as a source of energy, can be a useful indicator of potential changes in marine environments due to climate and ocean changes (Horwood et al., 2000;Hays et al., 2005).
In the face of such climate and ocean changes, we propose to develop EBRPs for exploited fish species and explore the ecosystem effects of fisheries exploitation under various scenarios of plankton productivity states (changes in plankton biomass) and fisheries management strategies with the latter consisting of two fishery scenarios (SS-focused and MS-focused), various fishing mortality rates, and two harvest polices (with and without harvest control rules, HCRs). We used an existing ecosystem model (OSMOSE-PNCIMA) for the Pacific North Coast Integrated Management Area (PNCIMA) off western Canada (Fu et al., 2017) to mimic the complex changes of climate and fisheries management strategies. The primary objectives of this study were to (i) derive ecosystem-level F MSY values for three exploited species [Pacific Herring (Clupea pallasii), Pacific Cod (Gadus macrocephalus), and Lingcod (Ophiodon elongatus)] under alternative scenarios of plankton productivity state and fisheries management strategy; (ii) explore the effects of different plankton productivity states and fishery management strategies on ecosystem-based F MSY ; and (iii) investigate ecosystem effects of interactions between plankton productivity states and fishery management strategies. The development of ecosystem-based BRPs in our study could be used to complement SS stock assessments by identifying key ecosystem dynamics that impact multiple species and fisheries over the long term, thus contributing to the need for incorporation of ecosystem considerations into scientific advice for sustainable fisheries management. Moreover, this study aims to guide the adaptation of fisheries management strategies in the face of climate and ocean changes.

Study area
The PNCIMA off western Canada encompasses around 88 000 km 2 , extending from the coastal watersheds to the outer limit of the continental slope. It 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;Fu et al., 2017). The spatial grid of the OSMOSE model developed in this study covers the whole PNCIMA region, divided into grid cells of 10 Â 10 km 2 ( Figure 1).

OSMOSE-PNCIMA model
The OSMOSE-PNCIMA model was originally developed by Fu et al. (2017), focusing on 6 key species, 19 "background" taxa, and 2 LTL groups (Supplementary Table S1). The key species are Pacific Herring, Pacific Cod, Lingcod, Arrowtooth Flounder (Atheresthes stomias), Walleye Pollock (Theragra chalcogramma), and Euphausiids (Thysanoessa spp. and Euphausia spp.). These key species are generally viewed as resident species among which strong predator-prey interactions occur (Fu et al., 2017). The inclusion of "background" taxa in OSMOSE-PNCIMA is to explicitly consider the species and taxa that are of secondary importance for this study but have potentially important trophic interactions with the modelled key species (see Supplementary  Table S1 for details of represented species in each background taxon). LTL groups (i.e. phytoplankton, and copepods) in OSMOSE-PNCIMA only serve as food and are represented as spatially distributed biomass pools (Supplementary Table S1; Fu et al., 2017).
Within the OSMOSE-PNCIMA model, Pacific Herring is assumed to have three separate populations in three distinct areas: the Prince Rupert District (PRD), the Haida Gwaii (HG), and the Central Coast (CC) based on genetic and tagging studies (Beacham et al., 2008;Flostrand et al., 2009). Pacific Cod and Lingcod have both been treated in stock assessments as having two different populations in two distinct areas: Hecate Strait (HS) and Queen Charlotte Sound (QCS) (King et al., 2012;Forrest et al., 2015a; Figure 1), albeit these distinctions are not based on genetic and tagging studies and simply use management areas as a proxy for stock structure. Therefore, a total of ten focus stocks were considered in OSMOSE-PNCIMA, i.e. Herring-PRD, Herring-HG, Herring-CC, Cod-HS, Cod-QCS, Lingcod-HS, Lingcod-QCS, Walleye Pollock, Arrowtooth Flounder, and Euphausiids. OSMOSE-PNCIMA simulated the life cycles of the ten focus stocks, from the egg stage to the terminal age, with a time-step of 3 months. At the first time-step after the production of eggs, the total number of eggs of each population was split into 120 super individuals called "schools," which were distributed spatially according to input distribution maps. The distribution maps (10 Â 10 km 2 ) were density-based and obtained from georeferenced data of both commercial fisheries and research surveys (data archives are maintained by Fisheries and Oceans Canada, Pacific Biological Station, Nanaimo, British Columbia). At each time-step, OSMOSE simulated the biological and ecological processes of these schools, including growth, predation, starvation, other natural mortality (due to causes unaccounted for by the model), fishing, reproduction, and spatial movement (including migration). However for the background taxa, only the predation, spatial distribution, and movement processes are simulated. Predation in OSMOSE was considered as an opportunistic process, occurring under the conditions of size suitability (with a minimum and a maximum predator to prey size ratio) and spatiotemporal co-occurrence between a predator and its prey (Fu et al., 2013(Fu et al., , 2017. Biomass estimates of Herring and Cod from the most recently available stock assessments (Forrest et al., 2015a;DFO, 2019) were used to validate the OSMOSE-PNCIMA ecosystem model because these two species have the most complete and reliable stock assessment results, and the model generated good fits to the assessed biomass time series for Herring and Cod. For detailed information on biological parameters, spatial distributions, and fishing as well as model construction and validation procedures, readers are referred to Fu et al. (2017).

Simulation scenarios
Changes in plankton biomass is often considered as a comprehensive indicator that accounts for climate and ocean changes in physical and biological functioning of marine ecosystems (Travers et al., 2007;Barange et al., 2014;Fu et al., 2018). Therefore, we used the change in plankton biomass (i.e. the sum of both phytoplankton and copepod biomass in the OSMOSE-PNCIMA model) as a proxy for plankton productivity change. Three plankton productivity states were thus hypothesized as current, high, and low plankton biomass scenarios. The assumed current plankton biomass scenario, denoted as Curr, was based Ecosystem-based reference points on the plankton biomass in the existing OSMOSE-PNCIMA model (Fu et al., 2017). The high plankton biomass (denoted as Doub) and low plankton biomass (denoted as Half) scenarios represented doubling and halving of the plankton biomass from the current level in the last 30 years, respectively. Due to the fisheries interaction effects (i.e. fishing one species may indirectly affect other species), which may affect the estimation of BRPs (Smith et al., 2015), this study explored two different fishery scenarios, single-species-focused (SS) and multi-species-focused (MS), as suggested by Smith et al. (2015). Nevertheless, BRPs estimated in both SS and MS scenarios were ecosystem-based BRPs which account for ecosystem processes. In the SS scenario, only one target species was subjected to fishing mortality experiments (i.e. a total of 33 fishing mortality rates changing from 0 to 2 year À1 were employed for each simulation run, fishing mortality rates were listed in Supplementary Table S2) while all other species were fished at the baseline fishing mortality rates. In the MS scenario, all the exploited key species undergo the same set of experimental fishing mortality rates from 0 to 2 year À1 .
Two harvest policies were applied to the focus stocks: HCR, denoting a policy with an applied HCR (details below), and NHCR, where no HCR was applied. For the HCR scenario, the fishing mortality rate of a specific species was subjected to a HCR based on the status of spawning stock biomass (SSB), following Fu and Schweigert (2004), i.e.
where t is for year, and SSB lowi and SSB up i are the lower and upper SSB cutoff and target thresholds for species i. Here, we used the same SSB thresholds for Pacific Herring, Pacific Cod, and Lingcod, i.e. SSB low ¼ 0.30*SSB 0 and SSB up ¼ 0.60*SSB 0 . The values of SSB low and SSB up used here were introduced for the potential management of Pacific Herring in British Columbia, Canada, and termed as limit reference point (LRP) and operational control points (OCPs), respectively . Here, LRP describes biomass level to be avoided with high probability; OCP specifies the biomass level where management action of reducing fishing mortality rate needs to take place . In total, 12 scenarios were simulated with the combination of three plankton productivity states (Doub, Curr, and Half), two fishery scenarios (SS and MS), and two harvest policies (HCR and NHCR) (see Supplementary Table S2).
Under each SS scenario, simulation runs were carried out using 1 of the 33 experimental fishing mortality rates ranging from 0 to 2 year À1 for the focus stocks (fish mortality levels were listed in Supplementary Table S2), while holding fishing mortality at the current levels for all other exploited species. Under the HCR scenario, the experimental fishing mortality was further subjected to the HCR according to the relative SSB [Equation (1)]. The experimental fishing mortality rates were applied to the focus stocks for the last 30 years of the simulation runs, after a 44-year burnin period, following Fu et al. (2017). The outputs of the last 10 years were summarized to calculate F MSY for all the focus stocks and to analyse the ecosystem effects (see following section). For each focus stock under the NHCR scenario, all the long-term catch values were plotted against experimental fishing mortality rates, and Locally Weighted Scatterplot Smoothing (Cleveland and Devlin, 1988) was fitted to the data points. Based on the constructed catch curves for each exploited stock (i.e. long-term fisheries catches as a function of fixed fishing mortality rates), F MSY was defined as the fishing mortality rate that achieved maximum yield (MSY) for the stock. This procedure of estimating F MSY is relatively standard and has been routinely employed in the field of ecosystem modelling (Smith et al., 2015;Grüss et al., 2016;Moffitt et al., 2016).

Indicators for assessing ecosystem effects
To investigate ecosystem dynamics under different scenarios, system catch (Y, represented by the sum of catches of all exploited species) and system biomass (B, represented by the sum of biomass for the same exploited species) were calculated and used as system-level indicators. In addition to the system catch and system biomass, the variability of these quantities can also function as indicators for ecosystem status and fisheries performance (Carpenter and Brock, 2006;Shin et al., 2010;Fu et al., 2013). The variability of B can be captured by its coefficient of variation (CV), denoted as CVB. This indicator reflects overall ecosystem stability (Shin et al., 2010) with a high CVB indicating a low ecosystem stability and low resistance to perturbations (Shin et al., 2010). The CV of system catch (CVY) was employed as an indicator of fisheries performance as achieving low inter-annual variability in fishery production is a desirable objective in fisheries management (Punt and Ralston, 2007;Punt, 2011). To improve comparability, both CVB and CVY were double square root transformed (Clarke and Warwick, 2001) and plotted against fishing mortality.

Estimation and comparison of ecosystem-based F MSY and MSY under various scenarios
Catch curves based on fisheries management without HCR (i.e. NHCR) yielded incrementally higher F MSY and MSY values for all seven focus stocks (i.e. Herring-PRD, Herring-HG, Herring-CC, Pacific Cod-QCS, Pacific Cod-HS, Lingcod-QCS, and Lingcod-HS) when plankton productivity shifted from low, current, to high scenario, respectively (Figures 2-8).
Under NHCR, it is worth noting that ecosystem-based F MSY values increased dramatically for all three herring stocks (i.e. Herring-PRD, Herring-HG, and Herring-CC) when fishery scenarios moved from SS-focused to MS-focused (Figures 2-4). Specifically, for Herring-PRD, the ecosystem-based F MSY values of 0.31, 0.40, and 0.42 in the low, current, and high plankton biomass scenarios, respectively, under the SS scenario (Figure 2     In terms of ecosystem-based MSY, the values estimated for all Pacific Herring (i.e. Herring-PRD, Herring-HG, and Herring-CC) and Pacific Cod (i.e. Pacific Cod-QCS and Pacific Cod-HS) stocks were relatively lower in SS than those in the MS scenario (Figures 2-6). However, MSY differences between SS and MS scenarios for two Lingcod stocks were only marginal with relatively higher MSY in the SS scenario (Figures 7  and 8).

Ecosystem effects of various plankton productivity states and fisheries management strategies
To investigate ecosystem effects under different plankton productivity states and fishing, we only focused on the MS fishery scenario. Comparisons of the ecosystem indicators in the SS fishery scenario were provided in Appendix (Supplementary Figures S1-S6).
As expected, both system catch and system biomass decreased as plankton productivity states changed from high, current, to low plankton biomass regardless of fisheries management with or without HCR (Figures 9 and 10). Specifically, long-term system catch curves under the NHCR scenario were dome-shaped under all plankton productivity states peaking at fishing mortality rates of 0.90, 0.75, and 0.55 in high, current, and low plankton biomass scenarios, respectively (Figure 9a, c, and e). However, under the HCR scenario, long-term system catch increased with the increase in fishing mortality in the high plankton biomass scenario ( Figure 9b). However, no obvious trend was seen in the current plankton biomass scenario (Figure 9d). In contrast, long-term system catch generally showed a decreasing trend when fishing mortality was >0.1 in the low plankton biomass scenario (Figure 9e). Although higher long-term system catch was achieved for fishing mortality rates ranging from 0.2 to 1 in the NHCR scenario compared to the HCR scenario, system catch in the HCR scenario was more consistent even under high fishing mortality >1.0 (Figure 9). In contrast with system catch, system biomass in the NHCR scenario was always lower than that in the HCR scenario ( Figure 10). System biomass in the NHCR scenario decreased with increasing fishing mortality under all three plankton productivity states (Figure 10a, c, and e). However, under the HCR scenario, system biomass remained stable even under high fishing mortality (Figure 10b, d, and f).
Both CVB and CVY increased as plankton productivity states varied from high, current, to low plankton biomass, regardless of whether HCR was applied (Figures 11 and 12). Under the NHCR scenario, the values of CVB and CVY generally increased with fishing mortality rates and the inter-stock variability (represented by the height of a box plot) of both CVB and CVY also increased with fishing mortality rates particularly under the high plankton biomass (Figure 11 NHCR and Figure 12 NHCR). In contrast, values of CVB and CVY in the HCR scenario were generally more stable along the gradient of fishing mortality. In addition, the inter-stock variability of both CVB and CVY were generally lower in the HCR scenario than under the NHCR scenario particularly under high fishing mortality >0.8 ( Figure

Discussion
With the end-to-end OSMOSE-PNCIMA model, this study estimated the ecosystem-based F MSY values for Pacific Herring, Pacific Cod, and Lingcod under three plankton productivity states. Main findings of this study include: (i) plankton productivity change (approximated by changes in plankton biomass) affected the values of ecosystem-based F MSY and MSY, which increased as plankton productivity states changed from low to high plankton biomass; (ii) ecosystem-based F MSY and MSY for Pacific Herring and Pacific Cod stocks increased when fishery scenarios shifted from SS to MS; and (iii) fisheries management incorporating HCRs helped maintain stable ecosystem production and sustainable fisheries.

Moving beyond SS-based BRPs
Developing BRPs is an essential component of decision-making for fisheries management (Restrepo et al., 1998;Gabriel and Mace, 1999;Sainsbury, 2008;Holsman et al., 2016). While the literature supports development of reference points that account for multispecies interactions (Gislason, 1999; Collie and Gislason  2001; Walters et al., 2005;Tyrrell et al., 2011;Holsman et al., 2016), in practice, reference points usually account only for SS population dynamics. However, advancement towards implementation of EBFM would benefit from development of EBRPs that account for key ecosystem processes and environmental changes (Link, 2010;Moffitt et al., 2016).
SS stock assessments were available for all the exploited stocks included in our OSMOSE-PNCIMA model, including estimates of BRPs (King et al., 2012;Forrest et al., 2015a;Grandin and Forrest, 2017;DFO, , 2019. For instance, F MSY was estimated at 0.298 year À1 and 0.312 year À1 for Pacific Cod-QCS and Pacific Cod-HS, respectively (higher than the ecosystem-based F MSY values estimated from OSMOSE-PNCIMA), although these were not adopted for management purposes (Forrest et al., 2015a). For Lingcod, F MSY was estimated at 0.07 year À1 and 0.09 year À1 for Lingcod-QCS and Lingcod-HS, respectively (lower than the ecosystem-based F MSY values estimated from OSMOSE-PNCIMA) (King et al., 2012). With the development of ecosystem models, it is now possible to explore ecosystem-based BRPs in fisheries management to better address issues associated with ecosystem processes and species interactions (Walters et al., 2005;Smith et al., 2015;Holsman et al., 2016).
Over the last decade, managers and fisheries scientists globally were committing to the development of methods for integrating ecosystem considerations and objectives into the estimation of ecosystem-based BRPs (Link, 2010;Fogarty and McCarthy, 2014;Buchheister et al., 2017). This study is the first attempt to investigate ecosystem-based F MSY for seven fish stocks in PNCIMA by using OSMOSE-PNCIMA (Fu et al., 2017). Although disparities exist between the previously estimated SS-based F MSY and our currently estimated ecosystem-based F MSY , such disparities were expected, and having also been noted by previous studies comparing SS and MS (or ecosystem-based) BRPs (Gislason, 1999;Holsman et al., 2016;Moffitt et al., 2016).

EBRPs under various plankton productivity and fishing scenarios
While BRPs were usually assumed to be equilibrium-based, there was growing recognition that BRPs for fisheries management were not static quantities, but instead may shift due to either internal or external factors (Murawski et al., 2001;Haltuch et al., 2009;Haltuch and Punt, 2011). For example, Collie and Gislason (2001) illustrated that life-history parameters as well as external factors, such as climate change, regime shift, eutrophication, or other environmental fluctuations could affect the values of BRPs. Chang et al. (2011) demonstrated that high-temperature-induced extra mortality under global warming could affect stock productivity and BRPs values, resulting in high risk of overexploitation in the long term if these changes were not incorporated into fisheries management. Holsman et al. (2016) suggested that the inclusion of environmental variables in projections, such as the variability in water temperature, may have a strong effect on the estimation of BRPs (see also in Basson, 1999). Moreover, Tyrrell et al. (2011) showed that incorporating predation mortality led to more conservative BRPs. A recent study using an EwE model built for Mille Lacs Lake also highlighted the trophic and temperature effects on the MSY reference points estimation (Kumar et al., 2017).
This study highlighted that ecosystem-based F MSY and MSY were affected by plankton productivity states. This was generally consistent with previous conclusions drawn by Heino et al. (2013) and Holsman et al. (2016), both of which showed the similar climatic effects on the estimation of BRPs. In addition, this study demonstrated that EBRPs were stock-specific in terms of the magnitude of difference in value among the three plankton productivity states. For instance, while EBRPs of Pacific Herring differed notably across the three plankton productivity states, those of Lingcod stocks did not show much difference among different plankton productivity states. Therefore, the impact of climate on EBRPs was stock-specific. Different variation of EBRPs in prey (e.g. Pacific Herring) and predator (e.g. Pacific Cod and Lingcod) stocks under various plankton productivity states may be due to the different trophic roles that these explored species play in the marine food web resulting in different trophic interactions (Travers-Trolet et al., 2014;Fu et al., 2017).
In addition, increasing values of F MSY and MSY for Pacific Herring and Pacific Cod stocks from SS to MS scenarios, may mainly be due to the effects of fisheries interactions which is consistent with the previous findings from an Atlantis model in Benguela system (Smith et al., 2015). In fact, exploitation of one species will invariably have some impact on other species in the system. Effective EBFM required an understanding of how a fishery that targets one species may indirectly affect other species in the same ecosystem. Competition and predator/prey interactions are important considerations in the design of sustainable fisheries management plans (Smith et al., 2015;Fu et al., 2017).

Effects of HCR in fisheries management
Fisheries management required the evaluation of harvest polices to evaluate trade-offs among multiple objectives such as meeting conservation concerns and achieving sustainable economic benefits (Punt, 2017). Harvest polices with cutoff thresholds (i.e. LRP) may be more robust to adverse environmental regimes under which the average survival rate falls below the long-term average for extended periods (Fu et al., 2000;King et al., 2014;Punt et al., 2014). To implement harvest polices with output controls (i.e. setting catch quotas), HCRs prescribe the annual fishing mortality as a function of stock status (Deroba and Bence, 2008). The purpose of implementing HCR was to reduce fishing mortality such that the stock biomass can grow to a target level where desirable economic outcomes are achieved and serious harm to the stock is avoided (Punt et al., 2008;Sainsbury, 2008;DFO, 2009). In short, HCRs were explicit guidelines designed to prevent future stock collapses, allow rebuilding depleted stocks, or maintain stocks at healthy levels (Kvamsdal et al., 2016).
Cutoff thresholds were often set at 25% of unfished equilibrium abundances or biomass (i.e. 0.25*SSB 0 , Thompson, 1993) and have been used for the management of Pacific Herring stocks in British Columbia, Canada since 1986 (Cleary et al., 2010). Alternatively, the Department of Fisheries and Oceans Canada's harvest policies apply provisional LRP and OCP of 0.4*B MSY and 0.8*B MSY , respectively (DFO, 2006(DFO, , 2009, and Kronlund et al. (2018) recently suggested a LRP of 0.30*SSB 0 for managing Herring Stocks in British Columbia, Canada. In an evaluation of common BRPs used in fisheries management, control rules based on estimated B 0 and stock depletion may perform better than those based on the estimates of B MSY , mainly because of the difficulty in estimating the latter (Haltuch et al., 2009;Forrest et al., 2018). Therefore in this study, we used 0.30*SSB 0 and 0.60*SSB 0 as LRP and OCP in the HCRs respectively (DFO, 2019), for

Ecosystem-based reference points
Pacific Herring, Pacific Cod, and Lingcod stocks. Results of our simulations generally showed that application of HCRs yielded higher system catch and system biomass as well as lower CVY and CVB, suggesting that fisheries management incorporating HCRs tended to maintain a more stable ecosystem and more sustainable fisheries.
However, different stocks responded quite differently to HCRs. For instance, HCRs were more effective for forage Pacific Herring stocks than for predatory Pacific Cod and Lingcod stocks. Different LRPs of HCRs have often been derived for different stocks depending on their productivity Kronlund et al., 2018), data availability (Forrest et al., 2015a;Holt et al., 2016), and sometimes environmental variability (Froese et al., 2011;Little et al., 2011). Thus, further investigation of system-and stock-specific HCRs for Pacific Cod, and especially Lingcod was warranted in future research. Moreover, evaluation of the performance of alternative HCRs for stocks in a feedback simulation context was also important for evaluating acceptable outcomes and trade-offs (Haltuch et al., 2009;Cleary et al., 2010;Forrest et al., 2018). More advanced simulations using the OSMOSE platform were underway for obtaining stock-specific thresholds of HCRs for these predatory stocks. Such ecosystem simulations can further explore alternative decisions on harvest rates and cutoff thresholds that explicitly consider the trade-offs among multiple ecosystem objectives, such as maintaining socioeconomic outcomes, biodiversity and ecosystem health (Fulton et al., 2014;Grüss et al., 2016).

Adapting fisheries management strategies in the face of climate and ocean changes
One of the core questions for fisheries scientists and managers was how to develop fisheries management strategies that are adaptive to uncertain environmental changes (Sumaila et al., 2011;King et al., 2014). For instance, Gaines et al. (2018) showed that improvement of fisheries management could offset many negative effects of climate change.
Our simulations compared the combined ecosystem effects of alternative plankton productivity states and fisheries management strategies, and may contribute to the development of adaptive fisheries management strategies in a changing environment for the PNCIMA ecosystem. For instance, when the plankton productivity state changed from low to high plankton biomass, the PNCIMA ecosystem tended to have higher system catch and system biomass but lower CVY and CVB, implying that more sustainable fisheries and stable ecosystem status was achieved. In addition, our simulation results generally concluded that to adapt to climate and ocean changes, incorporating suitable HCRs in fisheries management could benefit not only ecosystem health but also fisheries sustainability. The current study has taken a positive step forward by applying an ecosystem platform to facilitate the development of adaptive fisheries management strategies and the implementation of EBFM in a changing future.

Conclusions and future work
Although more progress is needed, ecosystem models can play a critical role in helping fisheries scientists to account for key ecosystem processes and environmental changes in fisheries management and provide guidance on adaptation to future changes in climate and human activities (Fulton, 2010;Fulton et al., 2014). With the application of the end-to-end OSMOSE model, this study has drawn the following conclusions that are potentially important for the future fisheries management: (i) ecosystembased BRPs complement SS BRPs by incorporating ecosystem considerations; and (ii) adopting HCRs can yield more stable ecosystem states and sustainable fisheries, especially in the face of climate and ocean changes.
This study, although not intended to provide tactical guidance regarding the fisheries management practices, helps build better knowledge of ecosystem-based BRPs and climate change impacts and may thus help identify future avenues towards EBFM. Future investigations may be aided by comparisons across a suite of models including minimum realistic models (Punt and Butterworth, 1995) and models of intermediate complexity (Plagányi et al., 2014) which represent attempts to introduce explicit interactions between multiple species while retaining the model precision necessary for tactical decision-making. Moreover, management strategy evaluation , a decision support tool that has been increasingly used in fisheries management for the last decade (Fulton et al., 2014;Grüss et al., 2016;Punt et al., 2016), can be conducted to determine the performance of the existing management strategy and to identify a "best" one among a set of candidate strategies. In addition, although plankton productivity employed in this study can be a useful forcing function for mimicking climate and ocean changes (Horwood et al., 2000;Hays et al., 2005), metabolic effects on fish from climate and ocean changes, such as water temperature, hypoxia, and ocean acidification, are also important considerations that have increasingly drawn attention over the last decade (Pörtner and Knust, 2007;Pörtner and Peck, 2010). Future simulations with ecosystem models may benefit from incorporating experimental results on metabolic impacts, directly or indirectly (e.g. reviewed in Koenigstein et al., 2016). Last but not least, this paper only considered a single climate driver, an extension of this study could evaluate synergistic and/or antagonistic effects of multiple drivers like ocean warming, ocean acidification and fishing (Griffith et al., 2012).

Supplementary data
Supplementary material is available at the ICESJMS online version of the manuscript.