Tradeoffs of managing cod as a sustainable resource in fluctuating environments

Sustainable human exploitation of living marine resources stems from a delicate 18 balance between yield stability and population persistence to achieve socioeconomic and 19 conservation goals. But our imperfect knowledge of how oceanic oscillations regulate temporal 20 variation in an exploited species can obscure the risk of missing management targets. We 21 illustrate how applying a management policy to suppress fluctuations in fishery yield in variable 22 environments (prey density and regional climate) can present unintended outcomes in harvested 23 predators and the sustainability of harvesting. Using Atlantic cod (Gadus morhua, an apex 24 predatory fish) in the Barents Sea as a case study we simulate age-structured population and 25 harvest dynamics through time-varying, density-dependent and -independent processes with 26 a stochastic, process-based model informed by 27-yr survey data. In this model, capelin 27 (Mallotus villosus, a pelagic forage fish), a primary prey of cod, fluctuations modulate the 28 strength of density-dependent regulation primarily through cannibalistic pressure on juvenile cod 29 survival; sea temperature fluctuations modulate thermal regulation of cod feeding, growth, 30 maturation, and reproduction. We first explore how capelin and temperature fluctuations filtered 31 through cod intrinsic dynamics modify catch stability and then evaluate how management to 32 suppress short-term variability in catch targets alters overharvest risk. Analyses reveal that 33 suppressing year-to-year catch variability impedes management responses to adjust fishing 34 pressure, which becomes progressively out of sync with variations in cod abundance. This 35 asynchrony becomes amplified in fluctuating environments, magnifying the amplitudes of both 36 fishing pressure and cod abundance and then intensifying the density-dependent regulation of 37 juvenile survival through cannibalism. Although these transient dynamics theoretically give 38 higher average catches, emergent, quasi-cyclic behaviors of the population would increase long39 term yield variability and elevate overharvest risk. Management strategies that overlook the 40


analysis
INTRODUCTION than three million t per yr (Dolgov 2002). Amplified year-to-year fluctuations in capelin in the 156 1980s and 1990s also shaped cod biomass rebuilding patterns: for example, cod cannibalism To simulate individual growth, the model converts the amount of consumed food (kcal) to 202 changes in body mass and length (cm) using temperature-and cod density-dependent functions 203 (Fig. 1d and Table 1: eq. 5a,b, Jones 1978). For simplicity we assume that the energy content of 204 consumed food is time-invariant and estimated through model fitting. The model then uses a 205 length-dependent function to simulate maturation for four and five-yr-olds and a length-and 206 condition-dependent function for six to nine-yr-olds ( Fig. 1d and Table 1: eq. 6a,b). The latter 207 function uses the amount of consumed food relative to growth as a proxy for condition (Fig. 1d 208 and Table 1: eq. 6c). All one to three-yr-olds are immature and all ten to fifteen-yr-olds mature 209 following the assumption adopted in the Arctic Fisheries Working Group stock assessment 210 (ICES 2011).

211
Parameterization-We used an integrated approach to parameterizing the submodels: 212 consumption, diet composition, growth, and maturation (Appendix S1: Fig. S1). This approach 213 uses reconstructed time series of age-specific numbers and fishing mortality rates (Appendix S1: 214 Table S1 and S2) taken from the stock assessment and survey data on prey biomass and sea 215 temperature (ICES 2011) as input to fit the model through the following four-step procedure 216 (Appendix S1: Fig. S1). First, we fit the model to total stomach content mass data to estimate 217 consumption model (Table 1: eq. 3a-e) parameters (β1-5, γ1-3, and ω, Table 2 and Appendix S1: 218 Table S5) using a nonlinear least squares (NLS) method (Appendix S1: Fig. S2). Second, we fit 219 the model to stomach content composition data to estimate diet composition model (Table 1: eq. 220 4a,b) parameters (DCothermin and ρ, Table 2) including age-specific suitability index for each prey 221 (S, Table 2 and Appendix S1: Table S6) using a correlation maximization method (Appendix S1: 222 Fig. S3). The first two steps were performed iteratively (Appendix S1: Fig. S1) until the model 223 captured the observed patterns. Third, we fit the model to cod mass data (Appendix S1: Table   224 S3) to estimate growth model (in mass and length, Table 1: eq. 5a,b) parameters (δ1-2 and μ1-3, 225 Table 2 and Appendix S1: Table S5) using NLS (Appendix S1: Fig. S4a). And forth, we fit the 226 model to maturity data (Appendix S1: Table S4) to estimate maturation model (Table 1: eq. 6a-c) 227 parameters (κ1-2 and ε1-3, Table 2 and Appendix S1: Table S5) using NLS (Appendix S1: Fig.  S4b). All tuning data used in this procedure had been collected during 1984-2004 (in which all 229 data were consistently available). 230 Cod management measures 231 We simulated dynamic cod mortality by fishing using a model based on the current harvest rule 232 set for Barents Sea cod (with some simplifications, see Appendix S1, ICES 2016) and evaluated is adjusted to Fpa scaled to the proportion of SSBy relative to Bpa (Table 1: eq. 7). Also, the 242 harvest rule has a policy tool to suppress between-year variability in Ctarget (yield stability 243 control) within ±10% (prior to the latest evaluation, ICES 2016). We computed a relative change 244 in Ctarget from year y-1 to y (interannual catch variability, ICVy, %, Table 1: eq. 8). In 245 simulations, when ICVy exceeds the set proportion, the algorithm searches for Ftarget,y until it 246 finds the largest value that does not exceed the bounds set by the stability control. But the 247 precautionary principle also applies to this policy. The stability control is applied only when 248 SSBy equals to or exceeds Bpa; this can result in asymmetrical applications (more frequent 249 applications in the positive direction). And the stability control is not applied in the first year 250 after the recovery period (SSBy ≤ Bpa). Once Ftarget,y is set, fishing mortality rates of a-yr-olds in 251 year y (Fa,y) are computed relative to those in the initial year of simulations (1990 in this study) 252 (  we randomly selected temperatures from these three groups sequentially; the length of each 264 period was randomly set so that the simulated thermal periods (cold, moderate, and warm) would 265 last from one to five years, effectively generating roughly decadal (autocorrelated) oscillations 266 (Fig. 1b).

267
Capelin production-We assume that capelin production (as change in total biomass) depends 268 on previous year's cod and capelin abundance to capture the historical patterns of density-269 dependent predation and autocorrelation (respectively). Although capelin is targeted for 270 commercial harvest in the system, the model does not simulate this fishery explicitly. We 271 empirically generated a variable production scenario using the historical  biomass 272 data (Appendix S1: Fig. S5) as described in Howell et al. (2013). Briefly, annual capelin biomass 273 was randomly selected with replacement from the dataset twice every year in simulations using 274 empirical functions: one based on cod SSB and another on capelin biomass in the previous year, 275 and then the mean was used as input (Howell et al .2013). To prevent biologically implausible 276 fluctuations in simulations, we fixed capelin biomass at six million t when projected cod SSB 277 exceeded 800,000 t. With this approach, simulated capelin production stochastically varied 278 among years and replicates (Appendix S1: Fig. S5), generating low-frequency oscillations (Fig. 279 1c).

280
Stochasticity in cod recruitment-To account for process uncertainty in cod productivity, we 281 added noise to simulated recruit numbers (Table 1: eq. 1). We randomly generated normally 282 distributed deviates, N(0, σ 2 ), where σ is the standard deviation estimated through model-fitting 283 using the SSB-recruit historical data as described above (see Food-web modeling); this process 284 was independently repeated every year for each replicate in simulations. Blim is set to 220,000 t (ICES 2016). 309 We initialized simulations with the 1990 estimates of age-specific demographic parameters and simulations with different replicate numbers remained < 0.06, Appendix S1: Fig. S6 and S7).

317
Because this study primarily focused on the implications of suppressing yield fluctuations, we 318 did not account for other sources of uncertainty such as population estimation and 319 implementation errors (we assume that annual yields are equal to Ctarget). We ran 325 simulations 320 in total to evaluate ecological assumptions (models) (25 x 4 = 100 simulations) and harvest 321 scenarios (25 x 9 = 225 simulations).

340
Ecological uncertainty in cod harvest dynamics 341 Density-dependent cannibalism played a key role in driving nonstationary patterns of cod 342 population and catch fluctuations under yield stability control in simulations (Fig. 2). Without 343 cannibalism, both varying capelin production and sea temperature (M1 and M2) induced low 344 variability in cod (adult and juvenile) biomass (less than 0.04 of the CV across years, Appendix 345 S2: Fig. 1a), and ICVs remained less than 2% (the stability control was never triggered, Fig. 2b).

346
These patterns reflect the extrinsic signals filtered through life-history processes such as 347 recruitment and maturation (but without intracohort interactions) in age-structured dynamics, as  and catch all increased; for example, catch was on average as much as 22% higher with a 10% 369 constraint ( Fig. 3a and Table 3). But because juveniles increased at lower rates than adults, 370 juvenile to adult biomass ratios declined by up to 14% with ICVs constrained ( Fig. 3a and Table   371 3). This trend reflects increasingly intensified density-dependent regulation of recruitment and  Also, the stability control would not be reapplied until the SSB rebuilds above Bpa, allowing 396 gradual increases in fishing pressure (Fig. 4a). This harvest measure thus can not only generate 397 asynchronies between the cycles of stock size and fishing pressure but also asymmetries in the 398 cyclicity of fishing pressure.

399
Stability-sustainability tradeoffs 400 Although average catches increased (by 10%-22%) with the strength of stability control ( Fig.   401 3a), overharvest risks also became increasingly greater owning primarily to high mean ICVs 402 (Fig. 5). When we set ICVs to 30% or higher, the probability of SSBs falling below Blim 403 remained less than 2% (Fig. 5). By contrast, when we constrained ICVs to less than 30%, SSBs 404 fell below Blim more often (up to 6.5%) because of greater amplitudes induced by the 405 unsynchronized stock-fishery cycles triggered by strict stability controls (Fig. 4).

407
Our results show that management to suppress short-term yield fluctuations without accounting 408 for the interplay between intrinsic and extrinsic dynamics can inadvertently intensify fishing 409 pressure, amplify long-term yield fluctuations, and destabilize stock-fishery dynamics, thereby    Table 1. List of the main equations used in the Barents Sea food-web model (StoCoBar) and cod management model. The description of variables and parameters of the equations are given in Table 2.

Description Equation
Population dynamics age 5-15 , +1 = −1, age 1-2 and 8-15  proportion of maximum food consumption by a-yr-old cod in year y FC max,a,y maximum food consumption by a-yr-old cod in year y B a,y a-yr-old cod biomass (kg) in year y λ a,y food availability index for a-yr-old cod in year y β 1, β 2, β 3, β 4, β 5 cod feeding model parameters estimated through model-fitting This study W a,y a-yr-old cod mass (g) in year y ICES (2011) for historical   ΔW a,y a-yr-old cod growth in mass (g) in year y FC kacl,a,y amount of consumed energy (kcal) by a-yr-old cod in year y, which is converted from FC a,y using conversion factors in Appendix S1: Table S7 L a,y a-yr-old cod length (cm) in year y IMR and PINRO databases δ 1 , δ 2 cod growth (in mass) model parameters estimated through model-fitting This study ΔL a,y a-yr-old cod growth in length (cm) in year y μ 1 , μ 2 , μ 3 cod growth (in length) model parameters estimated through model-fitting This study Fat a,y condition index for a-yr-old cod in year y FC' a,y consumed food relative to growth in mass for a-yr-old cod in year y κ 1 , κ 2 cod condition index parameters estimated through model-fitting This study D a,y proportion of adults in a-yr-old cod in year y ε 1 , ε 2 , ε 3 cod maturation model parameters estimated through model-fitting This study F 5-10,1 mean fishing mortality rate of five to ten-yr-olds (the dominant age classes in observed catches) in 1990 (0.274)

Baseline scenario simulations of cod population and harvest dynamics
Using the food-web model (configuration M4) under the baseline scenario (no stability control), cod population dynamics stabilized by 2030 with, on average, ~2.3 million t stock size, ~1.1 million t SSB, 694 million recruits, and 680,000 t catch (Appendix S2: Fig. S2), which all aligned closely with the historical population status and harvest of Barents Sea cod (ICES 2017). Overall, between-year relative changes in cod mass, and consumption rate increased with age, whereas those in maturity rate remained similar among maturing fish (four to nine-yr olds) (Appendix S2: Fig. S2 and S3). Cod diet, on average, comprised of 10% juvenile cod (more than 66% as one-yr-olds), 21% capelin, and 69% others (Appendix S2:    simulations. In local wavelet spectra, color gradient indicates red areas being higher power (intensity of periodicities) to blue areas being lower power, white areas indicate regions influenced by edge effects (outside the "cone of influence") in which inferences cannot be made, and y-axis is in the logarithms to the base 2. Artwork: Courtesy of the Integration and Application Network, University of Maryland Center for Environmental Science (ian.umces.edu/symbols/).   Figure S1. Wavelet power spectra computed for (a) annual catch forecasts, (b) total harvestable (three-yr-olds and older) biomass (stock size), (c) adult abundance (SSB), (d) juvenile survival (abundance of three-yr-olds), and (e) fishing mortality rate of Barents Sea cod projected with a stochastic, food-web model (StoCoBar configuration M4) under baseline, ±30%, ±20%, and ±10% catch constraints. Color gradients indicate red areas being higher power (intensity of periodicities) to blue areas being lower power, white areas indicate regions influenced by edge effects (outside the "cone of influence") in which inferences cannot be made, and y-axis is in the logarithms to the base 2. Figure S2. Wavelet coherence spectra computed for total harvestable biomass (stock size) and fishing mortality rate of Barents Sea cod projected with a stochastic, food-web model (StoCoBar configuration M4) under ±30%, ±20%, and ±10% catch constraints. Color gradients indicate red areas being higher power (intensity of periodicities) to blue areas being lower power, white areas indicate regions influenced by edge effects (outside the "cone of influence") in which inferences cannot be made, and y-axis is in the logarithms to the base 2. Arrows indicate the two series being in-phase (pointing right), the two series being anti-phase (left), stock size being leading (down), and fishing mortality rate being leading (up).