ICES Journal of Marine Science: Journal du Conseil

Edmunds, Peter, J;Doo, Steve, S;Carpenter, Robert, C

doi: 10.1093/icesjms/fsaa015pmid: N/A

Abstract In this study, fore reef coral communities were exposed to high pCO2 for a year to explore the relationship between net accretion (Gnet) and community structure (planar area growth). Coral reef communities simulating the fore reef at 17-m depth on Mo’orea, French Polynesia, were assembled in three outdoor flumes (each 500 l) that were maintained at ambient (396 µatm), 782 µatm, and 1434 µatm pCO2, supplied with seawater at 300 l h−1, and exposed to light simulating 17-m depth. The communities were constructed using corals from the fore reef, and the responses of massive Porites spp., Acropora spp., and Pocillopora verrucosa were assessed through monthly measurements of Gnet and planar area. High pCO2 depressed Gnet but did not affect colony area by taxon, although the areas of Acropora spp. and P. verrucosa summed to cause multivariate community structure to differ among treatments. These results suggest that skeletal plasticity modulates the effects of reduced Gnet at high pCO2 on planar growth, at least over a year. The low sensitivity of the planar growth of fore reef corals to the effects of ocean acidification (OA) on net calcification supports the counterintuitive conclusion that coral community structure may not be strongly affected by OA. Introduction Ocean acidification (OA), the consequence of the solution of rising atmospheric concentrations of CO2 in the oceans (Broecker, 1975), gained widespread attention as a global threat to marine ecosystems ∼2000 (Kleypas et al., 1999; Langdon et al., 2000). Since this realization, numerous experiments have been conducted in which phylogenetically diverse taxa have been incubated in tanks maintained at high seawater pCO2 (and reduced pH), thus revealing impaired biological functionality that most frequently has involved net calcification (Kroeker et al., 2010). For calcified organisms, net accretion typically is reduced by lowered seawater pH (Chan and Connolly, 2012; Kroeker et al., 2013; Kornder et al., 2018), leading to the conclusion that ecosystems dominated by foundation species having calcified shells and skeletons will be threatened by OA (Doney et al. 2009). Tropical coral reefs conform to this prediction (Doney et al., 2009; Kornder et al., 2018), with the most severe of the possibilities (i.e. net ecosystem dissolution) predicted to occur within the current century on sediment-dominated reefs (Eyre et al., 2018). Much is known of the short-term effects of elevated seawater pCO2 on the net calcification of corals and coral reef communities (Chan and Connolly, 2012; Kroeker et al., 2013; Kornder et al., 2018). However, less is known about the response of coral reef communities to high pCO2 over ecologically relevant time scales of a year or more, with only a few studies conducted in the laboratory (e.g. Horwitz et al., 2017; Edmunds et al., 2019a) or the field (Fabricius et al., 2011). The limited number of studies addressing the response of coral communities to OA over long time scales creates a significant knowledge gap in understanding the effects of OA on coral reefs. Addressing this gap requires information on aspects of the biological effects of OA that have yet to be addressed in detail. One of these effects involves the processes by which net calcification influences the planar spread of scleractinian skeletons across the benthos, which represents the state variable (i.e. coral cover) through which coral community ecology has advanced for decades (Loya, 1972; Connell, 1973; De’ath et al., 2012). Indeed, the coral reef crisis largely is defined by the decline in coral cover (Bellwood et al., 2004; Hughes et al., 2010) and the gravity of future projections is evaluated by how little coral cover will remain within a century (Edmunds, 2015; Donner et al., 2018), together with low, or negative, rates of net community calcification (van Hooidonk et al., 2014; Eyre et al., 2018). There is a great need, therefore, to express the implications for coral reefs of OA using the same parameter (i.e. coral cover) through which their ecological success typically is evaluated. Testing for the effects of OA on the benthic ecology of coral reefs can only be accomplished through analyses of natural communities exposed to elevated seawater pCO2 (e.g. volcanic seeps, Fabricius et al., 2011), free ocean CO2 enrichment experiments (FOCE) (Kline et al., 2012; Doo et al., 2019), or laboratory experiments using tanks with the volumetric capacity to accommodate reef communities (Langdon et al., 2000; Jokiel et al., 2008; Carpenter et al., 2018). Analyses of coral reefs naturally enriched by CO2 are subject to the limitations of mensurative experimentation (Hurlbert, 1984; Cornwall and Hurd, 2015), while the ability of FOCE and laboratory experiments to test for effects using ecological response variables (e.g. planar area) requires incubations lasting long enough to allow slow-onset effects to be detected. Despite these challenges, several studies have begun to suggest that negative effects of OA on the ecology of coral reef communities are complex and subtle, even while their net calcification declines (Edmunds et al., 2019a; Rippe et al., 2018). Dove et al. (2013), for example incubated coral communities from Heron Island, Australia, for ∼84 d under four combinations of pCO2 and temperature using 300-l tanks. At 572 µatm pCO2 and +4°C, net community calcification shifted to dissolution and coral cover declined from ∼35% to ∼27%. Horwitz et al. (2017) conducted a year-long experiment in 30-l tanks and compared the effects of 400 vs. 1795 ppm CO2 on growth and competition among corals from Eilat, Israel. Coral–coral competition reduced growth to such an extent that there was no further additive effect of high pCO2 on coral growth. High pCO2 altered the competitive hierarchy among six coral species, which was inferred to promote shifts in community structure (Horwitz et al., 2017). Close spacing among coral colonies, also has the potential to modulate the effects of OA on coral growth, as aggregates of Pocillopora verrucosa are less affected by ∼1200 µatm pCO2 than single branches (Evensen and Edmunds, 2017). Moreover, multi-branched coral colonies are less susceptible to high pCO2 than single branches (Edmunds and Burgess, 2016). Reef communities in Papua New Guinea naturally exposed to ∼750 ppm pCO2 show changes (but not complete losses) in scleractinian communities relative to reefs at ∼390 ppm pCO2 (Fabricius et al., 2011), and in Japan, reefs are dominated (∼50% cover) by octocorals at 830 µatm pCO2 but depleted of octocorals and scleractinians (∼5% cover) at 1465 µatm pCO2 (Inoue et al., 2013). Working with back reef communities in Mo’orea, French Polynesia, Comeau et al. (2015) showed that net community calcification was reduced by 59% over 8 weeks at 1300 µatm pCO2, mostly through sediment dissolution. However, when an experiment with back reef communities that was comparable to that of Comeau et al. (2015) was extended over 1 year (employing pCO2 as high as 1067 µatm), this effect did not translate into changes in community structure, even though the net calcification of individual colonies was depressed (Edmunds et al., 2019a). By revealing the effects of OA on corals that emerge as the scale of investigation changes, for example from coral fragments to coral communities, the aforementioned studies highlight the challenges of inferring ecological consequences from the negative effects of OA on coral calcification and coral reef biogeochemistry. Responding to these challenges provided motivation for extending our analysis of the effects of OA on back reef community structure (Edmunds et al., 2019a) to a coral reef habitat (i.e. the fore reef) that is more complex in terms of species assemblages and physical environmental conditions. The fore reef habitat is vastly different from the back reef in terms of community assemblages (e.g. coral cover and diversity, Done, 1983) and physical environmental conditions (e.g. wave forces, light, and often temperature, Done, 1983). The large aerial extent of the fore reef habitat relative to back reef habitat, and the high rates of net accretion achieved by fore reef communities (Vecsei 2001, 2004), makes them a critical focus for studies of the effects of OA. Unlike our early year-long analysis (Edmunds et al., 2019a), here, we focus on fore reef communities from 17-m depth and extend our 7-week analysis (Comeau et al., 2016) of the effects of high pCO2 on fore reef communities to a year. The experiment tests the null hypothesis that community structure is unaffected by long-term exposure to elevated pCO2 (targeted at 700 and 1300 µatm). By simultaneously measuring net organismic calcification (Gnet), and community structure through two-dimensional planar images, we explored the relationship between a common measure of the effects of OA on coral reef organisms and a common community-based metric by which the ecological consequences of the coral reef crisis are expressed. Material and methods Overview Fore reef communities were assembled in three outdoor flumes in Mo’orea, which were assigned randomly to pCO2 treatments targeting ambient (400 μatm) and 700 and 1300 μatm pCO2. The elevated pCO2 treatments approximated atmospheric conditions projected for ∼2140 under representative concentration pathways 2.6, 4.5, and 8.5 (IPCC, 2014), respectively. Treatments were maintained for 1 year beginning in late Austral spring (November 2017), and actual pCO2 treatments over the year differed from target values (described below). The flumes are described elsewhere (Carpenter et al., 2018, Edmunds et al., 2019a). In brief, each flume consisted of a working section that was 5.0-m long and 30-cm wide and filled to ∼30-cm depth with ∼500 l of seawater. The working section contained fore reef community members that either were secured on the floor of the flumes (the “fixed” portion) or were placed on the floor of the flumes (i.e. they were “unfixed”), and together these portions occupied 4.7 m × 0.3 m of the floor of the working section of each flume. Seawater was circulated continually through a return section and was exchanged with fresh seawater at ∼300 l h−1. Seawater was pumped from Cook’s Bay (14-m depth) and filtered through sand (pore size ∼450–550 µm) before entering the flumes. With this pore size, small particulates passed through the filter and were added to the flumes where they were available as food for heterotrophic organisms. It was not logistically feasible to pump seawater from the fore reef (1.5 km away), but we reasoned that the high flow rate of fore reef water over the reef crest (Hench et al., 2008), and the pumping of water from 14-m depth would reduce artefacts related to the source of seawater used to fill the flumes. Our experiment is pseudoreplicated (both “simple” and “temporal”, after Hurlbert 1984), but we sought to balance the value of statistical independence against ecological relevance in terms of volume, functional scale, environmental conditions, and experimental duration. The challenges of achieving ecological relevance cannot always be reconciled with statistical perfection (e.g. Kline et al., 2012; Albright et al., 2018), but the experiment can be implemented to alleviate the lack of independence of replication that is the core concept of pseudoreplication. In the present study, the small size of the corals relative to the volume of seawater in the flumes, and the high exchange rate of seawater (∼60% h−1), reduced the likelihood that corals in any one flume were influenced by other corals in the same flume (i.e. they probably were independent). The inferred independence of corals within each flume reduced the implications of pseudoreplication in the analysis of coral growth (change in the area) over time, in which corals were treated as replicates. In the analysis of the net change in mass of corals, the dependence of measurements on single corals at multiple times was addressed with a repeated measures statistical design (both described in “Statistical Analyses” section). Fore reef communities The reef communities were assembled to correspond to the mean percent cover of the major benthic space holders recorded in 2006 at 17-m depth on the fore reef of the north shore of Mo’orea (Carpenter, 2018; Edmunds, 2019). A historic community structure (rather than present-day) was used because 2006 represented the long-term community structure on this reef (Edmunds et al., 2019b), and it created the capacity to compare aspects of the present experiment with a previous experiment (Comeau et al., 2016). Based on six sites sampled around Mo’orea in 2006, the community structure in the flumes was targeted to ∼11% cover of Pocillopora spp., ∼8% massive Porites spp., 8% Acropora spp., and ∼53% reef rock. This construct created a community with ∼27% coral cover, which was slightly lower than the actual mean coral cover in 2006 (32%) because the remaining 14 genera of scleractinians and Millepora contributed 5% coral cover. The Pocillopora conformed to the classic morphology of P. verrucosa (Veron, 2000), but it is likely that other Pocillopora spp. were present in the flumes (Edmunds et al., 2016). Likewise, Acropora spp. were selected to represent Acropora hyacinthus and A. retusa, which were common on the fore reef when the experiment was completed, and colonies of these species were scattered haphazardly among the flumes. Given the morphological complexity of Acropora spp., it is possible that other species were placed into the flumes. Pieces of coral rubble (∼11.5-cm diameter) were added to achieve ∼29% cover. Coral and rubble were haphazardly scattered along the working section of each flume to approach the targets for percentage cover, and this resulted in portions of the flumes having slightly different covers of coral. This was important for the central 2.4-m portion of the flume, where community members were fixed to allow the community structure to be quantified monthly using planar photographs. In the adjacent portions of the flumes, community members were unfixed so that they could be removed monthly for buoyant weighing (described below). Corals and rubble were collected from ∼17-m depth on the north shore fore reef, epoxied (Z-Spar A788; Pettit Marine Paint, Rockaway, NJ, USA) to plastic bases, and placed in a seawater table for at least 2 d before being added to the flumes. This time allowed the epoxy to cure and for the corals to recover from the collection. Fore reef communities were assembled in the flumes on 27 October 2017, where they were maintained under ambient seawater conditions until 3 November. At this time, treatment pCO2 levels were initiated in two flumes (one remained at ambient pCO2), with pCO2 gradually increased to target values over 24 h. Physical and chemical parameters Seawater was circulated in the flumes at ∼0.1 m s−1 using a pump (Wave II 373 J s−1, W. Lim Corp., El Monte, CA, USA), and flow speeds were measured across the working sections using a Nortek Vectrino Acoustic Doppler Velocimeter. This flow speed was ecologically relevant for 15-m depth on the fore reef of Mo’orea (14-year mean = 0.065 m s−1, Washburn, 2018). The flumes were exposed to natural sunlight that was reduced with a blue filter (LEE # 183; Lee Filters, Andover, England) to photon flux densities (PFD) in the range of photosynthetically active radiation (400–700 nm) that approximated those at 17-m depth (Comeau et al., 2016). Light in the flumes was measured continuously (at 0.0006 Hz) using cosine-corrected sensors (Odyssey; Dataflow Systems Ltd, Christchurch, New Zealand) that recorded photosynthetically active radiation (PAR). Odyssey sensors were calibrated with a Li-COR meter [LI-1400, Li-COR Biosciences, Lincoln, NE, USA attached to a 2π sensor (LI 192A)]. Temperatures in the flumes were regulated with chillers (heaters were not required) and were maintained close to the mean monthly seawater temperature at 17-m depth on the fore reef. Seawater carbonate chemistry was uncontrolled in one flume (ambient, ∼400 μatm pCO2), and controlled in two others to simulate conditions arising from seawater pCO2 targeted at 700 and 1300 μatm. Seawater pH was not altered in the ambient flume but was controlled in the treatment flumes by bubbling CO2 into the seawater to alter pH relative to a set point (regulated using an AquaController, Neptune Systems, Morgan Hill, CA, USA) that operated a solenoid supplying pure CO2 gas to a diffuser stone submerged in each flume. A diurnal downward pH adjustment of ∼0.1 unit was applied to the two treatment flumes to simulate natural diurnal variation in seawater pCO2 on the reef of Mo’orea (Hofmann et al., 2011). The ambient flume also maintained a diurnal variation in pCO2 with a night-time pH ∼0.1 unit lower than in the daytime. Ambient air was bubbled continuously into all flumes. Periodic measurements of pCO2 in the flumes confirmed that nocturnal pCO2 met, or exceeded, day-time target values (described in results). Throughout the experiment, logging sensors (described above) recorded PAR, and temperature [Hobo Pro v2 (±0.2°C); Onset Computer Corp., Bourne, MA, USA]. pH was measured daily on the total hydrogen ion scale (pHT) using a handheld meter (see below). The values from the temperature and pH measurements were used to adjust the thermostat and pH set points to achieve target pCO2 values. Seawater carbonate chemistry was calculated weekly using pH and AT measurements recorded once during the day, and once at night. A bench-top conductivity meter (Orion Star A212; Thermo Scientific, Waltham, MA, USA) was used to measure the salinity of the same water samples. The parameters of the seawater carbonate system were calculated from temperature, salinity, pHT, and AT, using the R package Seacarb (Lavigne and Gattuso, 2013). Calculations were made using the carbonic acid dissociation constants of Lueker et al. (2000), the KSO4 concentration for the bisulphate ion from Dickson et al. (1990) and the Kf constant of Perez and Fraga (1987). pHT was measured using a DG 115-SC electrode (Mettler Toledo, Columbus, OH, USA) that was calibrated with a TRIS buffer (SOP 6a, Dickson et al., 2007). AT was measured using open-cell, acidimetric titration (SOP 3b, Dickson et al., 2007) using certified titrant with an automatic titrator (T50, Mettler Toledo) fitted with a DG 115-SC electrode (Mettler Toledo). The accuracy and precision of measurements were determined by processing certified reference materials (CRMs’ batch numbers 158 and 172; A. Dickson Laboratory, Scripps Institution of Oceanography, CA, USA), against which measured values of AT maintained an accuracy of 2.7 ± 0.4 μmol kg−1 (n = 54) and the precision of 1.8 ± 0.1 μmol kg−1 (n = 475). Response variables Net changes in mass (Gnet) of corals in the unfixed portion of the community were measured every month by buoyant weighing (accuracy ±1 mg CaCO3) (Spencer Davies, 1989). The fixed community members were weighed at the start and end of the experiment. Buoyant weight was converted to dry weight of CaCO3 using empirical seawater density [1.02278 ± 0.00475 g cm−3 (mean ± s.e. based on 17 determinations over the year)] and the density of pure aragonite (2.93 g cm−3). As the three-dimensional area of tissue changed throughout the year as a result of growth and partial mortality, the change in mass could not be normalized to the actual surface area of the live coral tissue. Gnet at each time, therefore, was expressed as the percentage change in mass relative to the initial mass in November 2017. Community structure The effects of the treatments on the community structure were described using photographs recorded monthly in the planar view. The image-based technique strengthened the ability to address the effects of OA on the community ecology of coral reefs, which frequently is recorded using planar photographs (including in Mo’orea, Edmunds, 2019). Given the high capacity for corals on the fore reef of Mo’orea to increase in planar cover within a year (e.g. Edmunds, 2019), and the linear extension rates that can be expected for the corals employed in our study (e.g. https://coraltraits.org/), it was reasonable to expect changes in planar area of corals, and the effects of treatments on these areas, to be detected by photography over a year (e.g. Edmunds et al., 2019a). Photographs were recorded in ambient light using a GoPro Hero 4 camera (12 MP, 3-mm focal length) that was fitted to a stand and positioned on the upper edge of the flumes to record the benthic community through the air–water interface. At each sampling, the camera was sequentially moved along the flume to record the community in the middle 2.4 m of the working section using ∼15 contiguous photographs. Photographs were analysed using ImageJ software (Abramoff et al., 2004) after they were stitched together to make a single image for each sampling. This image covered the ∼2.4-m length of the central portion of the flume where the corals were secured to a plastic-coated metal grid with a mesh size of 5 cm × 5 cm (see Edmunds et al., 2019a). The stitching of photographs sometimes was imperfect due to parallax errors, and in such cases, separate pictures were evaluated to assess organism size. The planar area of living tissue on corals was quantified by outlining organisms in ImageJ, after scaling the image using the metal grid as a size reference. Organism size (cm2) was expressed as a percentage of the area (7200 cm2) occupied by the fixed portion of the community. The summed area of community members was used to determine the cover of the benthic community, and the areas of each organism were used to quantify growth (and shrinkage). Where organisms died, their area was set to zero. Statistical analyses Mortality of the corals was compared between the fixed and unfixed portions of the flumes, and among treatments, using χ2 contingency tables. Gnet was analysed by species using repeated measure permutational multivariate analysis of variance (PERMANOVA) (999 permutations, α = 0.05) in which the percentage change relative to the initial mass was the response variable, flume was the fixed effects, and time was the repeated measure. Data were square-root transformed and prepared as a resemblance matrix using Bray–Curtis dissimilarity values. The rate of change in the planar area for each organism was evaluated by the least squares linear regression of non-zero areas (i.e. when the organisms were alive) on time (months). Growth rates (i.e. the slopes of the regression lines, % month−1) were used to compare the planar growth of each taxon among flumes using Kruskal–Wallis tests. Multivariate community structure over time was described with two-dimensional ordinations using non-metric multidimensional scaling (nMDS). The resemblance matrices used for nMDS cannot accommodate missing data, and missing data either were set to zero (where organisms died) or replaced by values interpolated by the least squares linear regressions prepared for each organism using all available data for that organism. Resemblance matrices were based on percentage cover and were prepared independently for each flume. Data were log(x + 1) transformed and converted to Bray–Curtis dissimilarities; nMDS plots were obtained using 100 iterations until stress was <1.0. To test for differences among pCO2 treatments in the temporal variation in multivariate community structure, resemblance matrices were compared using rank correlation coefficients in a pairwise fashion against the resemblance matrix for the community incubated under ambient pCO2; significance was evaluated in a permutational framework (999 iterations) using an α of 0.05. The roles of chemical and physical conditions (Supplementary Table S1) in driving changing community structure were evaluated using multivariate tests of association for each flume. The BEST routine (in PRIMER, Clark and Gorley, 2006) was used to test for associations between community structure (evaluated using Bray–Curtis dissimilarities) and the chemical and physical conditions (evaluated using Euclidian Distances). Chemical and physical data (Supplementary Table S1) were screened for co-linearity, and one member of each pair of co-linear variables was excluded from the analysis. Univariate statistics were completed using SYSTAT 13 (SYSTAT Software Inc., San Jose, CA, USA), and multivariate and permutational statistics were completed using Primer 6.0 with the PERMANOVA+ add on (Clark and Gorley, 2006). Results Overview The experiment ran for 12 months over which the mean monthly temperature (MMT) ranged from 27.9 ± 0.2°C (1434 µatm flume) to 28.1 ± 0.2°C in the other treatments (± s.e., n = 12 months) and MMTs varied from 28.7 ± 0.1°C (January, March, and April) to 26.9 ± 0.1°C October (± s.e., n = 3 flumes). PFD in each flume was affected by daily variation in weather, and maximum daily PFD varied from 14 to 1636 μmol quanta m−2 s−1, with the highest values occurring when the lids of the flumes were removed during cleaning. Overall, mean (± s.e., n = 328 d) maximum irradiances were 370 ± 12 μmol quanta m−2 s−1 (400 µatm flume), 296 ± 7 μmol quanta m−2 s−1 (700 µatm flume), and 274 ± 9 μmol quanta m−2 s−1 (1300 µatm flume). The mean monthly maximum daily PFD ranged from 63 ± 8 µmol quanta m−2 s−1 (February 2018) to 641 ± 145 µmol quanta m−2 s−1 (June 2018), and daily integrated values ranged from 1.1 ± 0.1 mol quanta m−2 d−1 (February 2018) to 14.5 ± 1.9 mol quanta m−2 d−1 (September 2018). The three flumes contrasted mean day-time pCO2 treatments of 375 ± 15 (control), 711 ± 15, and 1315 ± 20 µatm, which corresponded to mean pHTotal and AT of 8.1 ± <0.1 and 2340 ± 2 µmol kg−1, 7.8 ± <0.1 and 2344 ± 3 µmol kg−1, and 7.6 ± <0.1 and 2349 ± 3 µmol kg−1, respectively (all ± s.e., n = 47, Supplementary Table S2). The night-time conditions of the three flumes were 416 ± 8 (control), 853 ± 23, and 1553 ± 30 µatm, which corresponded to mean pHTotal and AT of 8.03 ± <0.1 and 2336 ± 3 µmol kg−1, 7.76 ± <0.1 and 2338 ± 4 µmol kg−1, and 7.54 ± <0.1 and 2350 ± 3 µmol kg−1, respectively (all ± s.e., n = 40; Supplementary Table S2). Salinity values during the incubation were 35.4 ± <0.1 PSU, 35.5 ± <0.1 PSU, and 35.6 ± <0.1 PSU for the 400, 800, and 1300 µatm flumes, respectively (mean ± s.e., n = 47), and did not differ significantly from day to night (Supplementary Table S2). In November 2017, the working section of each flume contained 9–10 colonies of Acropora spp., 14–17 colonies of massive Porites spp., and 14–17 colonies of P. verrucosa, of which 4–5, 7–8, and 7–8 colonies, respectively, were in the fixed (2.4-m long) community section. Initial coral cover was 20% (in the 396 µatm flume), 19% (782 µatm), and 20% (1434 µatm) in the working section of each flume, with 18%, 20%, and 18% in the fixed community sections. This cover represented 4–7% of Acropora spp., 5–7% of massive Porites spp., and 6–9% of P. verrucosa. One year later, mortality left 9 colonies of Acropora spp., 1–3 colonies of massive Porites spp., and 11–15 colonies of P. verrucosa in each flume, with 4–5, 0–2, and 5–6 colonies, respectively, in the fixed community sections of each flume. The coral cover at the end of the experiment (November 2018) was 14% (396 µatm), 13% (782 µatm), and 12% (1434 µatm), with the fixed communities having a combined coral cover of 11%, 13%, and 11%, respectively. Net changes in mass (Gnet) In November 2017, the mean (± s.e.) mass of the corals in the unfixed portions of the communities was 77 ± 6 g (n = 25) for P. verrucosa, 104 ± 12 g for Acropora spp. (n = 15), and 95 ± 10 g for massive Porites spp. (n = 24); the fixed corals had similar masses [86 ± 7 g (n = 23), 116 ± 13 g (n = 13), and 86 ± 8 g (n = 23), respectively]. Over the year, a few Acropora spp. died (7%, n = 25 across all treatments), but mortality was higher for P. verrucosa (14%, n = 49) and highest for massive Porites (87%, n = 47). Mortality was independent of treatment for massive Porites spp. (χ2 = 1.267, df = 2, p = 0.531) and P. verrucosa (χ2 = 0.572, df = 2, p = 0.768) and independent of whether they were in the fixed or unfixed portion of the communities (pooled among taxa and treatments, χ2 = 0.654, df = 1, p = 0.418). By November 2018, among the unfixed members of the community, 96% of the corals in the 396 µatm treatment increased in mass (all of the Acropora spp. and massive Porites and 89% of the P. verrucosa), 74% of the corals in the 782 µatm treatment increased in mass (all of the P. verrucosa, 80% of the Acropora spp., and 50% of the massive Porites), and 68% of corals in the 1434 µatm treatment increased in mass (80% of the Acropora spp., 90% of the P. verrucosa, and 29% of the massive Porites). Among the fixed members of the community, 90% of the corals in the 396 µatm treatment increased in mass (all the Acropora spp. and massive Porites and 75% of the P. verrucosa), 91% of the corals in the 782 µatm treatment increased in mass (all of Acropora spp. and P. verrucosa and 75% of the massive Porites), and 83% of the corals in the 1434 µatm treatment increased in mass (all of the Acropora spp., 86% of the P. verrucosa, and 71% of the massive Porites). Gnet of the unfixed corals that were weighed monthly varied among taxa and treatments (Figure 1). To statistically test for effects of treatments and time, PERMANOVA requires complete data frames and missing data, therefore, were replaced with values interpolated by linear regression, with 4 of 195 values interpolated for Acropora spp., 7 of 312 for massive Porites spp., and 8 of 312 for P. verrucosa. For all three taxa, Gnet varied over time and among flumes, but there was no interaction between these effects (Table 1, Figure 1). Figure 1. Open in new tabDownload slide Summary of changes in net mass (Gnet) of corals incubated in flumes maintained at three different pCO2 treatments. Values show members of the unfixed community (open boxes, monthly determinations) and fixed community (shaded, November 2018 only) and display dry mass as a percentage of the initial mass of each coral. The percentage scale supports contrasts among corals that differed in initial mass. Box plots report medians as a line within each box (linked by red lines for the unfixed community), boxes show quartiles, and whiskers display ×1.5 the interquartile range, with outliners plotted individually. Sample sizes are shown in each plot frame for unfixed and fixed organisms. Figure 1. Open in new tabDownload slide Summary of changes in net mass (Gnet) of corals incubated in flumes maintained at three different pCO2 treatments. Values show members of the unfixed community (open boxes, monthly determinations) and fixed community (shaded, November 2018 only) and display dry mass as a percentage of the initial mass of each coral. The percentage scale supports contrasts among corals that differed in initial mass. Box plots report medians as a line within each box (linked by red lines for the unfixed community), boxes show quartiles, and whiskers display ×1.5 the interquartile range, with outliners plotted individually. Sample sizes are shown in each plot frame for unfixed and fixed organisms. Table 1. Summary of statistical contrasts of Gnet among treatments (i.e. flumes) and times for corals incubated for a year in three flumes maintained at different pCO2 levels. Taxon . Effects . Pseudo-F . df . p-perm . Acropora spp. Time 4.504 12 180 0.001 Treatment 40.676 2 180 0.001 Massive Porites spp. Time 1.876 12 297 0.031 Flumes 29.414 2 297 0.001 Pocillopora verrucosa Time 10.565 12 297 0.001 Flumes 19.968 2 297 0.001 Taxon . Effects . Pseudo-F . df . p-perm . Acropora spp. Time 4.504 12 180 0.001 Treatment 40.676 2 180 0.001 Massive Porites spp. Time 1.876 12 297 0.031 Flumes 29.414 2 297 0.001 Pocillopora verrucosa Time 10.565 12 297 0.001 Flumes 19.968 2 297 0.001 Analyses completed using repeated measures PERMANOVA in which flume was the fixed effects (3 levels), time was the repeated measures effects (13 levels), and relative Gnet was the dependent variable (and was square-root transformed data). Open in new tab Table 1. Summary of statistical contrasts of Gnet among treatments (i.e. flumes) and times for corals incubated for a year in three flumes maintained at different pCO2 levels. Taxon . Effects . Pseudo-F . df . p-perm . Acropora spp. Time 4.504 12 180 0.001 Treatment 40.676 2 180 0.001 Massive Porites spp. Time 1.876 12 297 0.031 Flumes 29.414 2 297 0.001 Pocillopora verrucosa Time 10.565 12 297 0.001 Flumes 19.968 2 297 0.001 Taxon . Effects . Pseudo-F . df . p-perm . Acropora spp. Time 4.504 12 180 0.001 Treatment 40.676 2 180 0.001 Massive Porites spp. Time 1.876 12 297 0.031 Flumes 29.414 2 297 0.001 Pocillopora verrucosa Time 10.565 12 297 0.001 Flumes 19.968 2 297 0.001 Analyses completed using repeated measures PERMANOVA in which flume was the fixed effects (3 levels), time was the repeated measures effects (13 levels), and relative Gnet was the dependent variable (and was square-root transformed data). Open in new tab Community structure At the start of the experiment, the size of the corals in the fixed community varied from 16 to 179 cm2. Most Acropora spp. grew, most P. verrucosa shrank, and most massive Porites spp. shrank and died (Figure 2, Supplementary Figure S1); mortality of massive Porites spp. rapidly began in month 5 of the experiment (March 2018) when seawater temperature was reaching its maximum annual value (Supplementary Table S1). Based on the regressions of size on time (months) (Supplementary Table S3), the rates of change in size (i.e. the slopes of these regressions) were unaffected by treatments for Acropora spp. (H = 0.119, n1 = 4, n2 = 5, n3 = 4, p = 0.942), P. verrucosa (H = 1.068, n1 = 7, n2 = 8, n3 = 9, p = 0.396), and massive Porites spp. (H = 1.851, n1 = 7, n2 = 7, n3 = 8, p = 0.396). Figure 2. Open in new tabDownload slide Mean (± s.e. confidence belts) percentage cover of live tissue on colonies of Acropora spp., Pocillopora verrucosa, and massive Porites spp. in the fixed portions of three flumes incubated at different pCO2 treatments (A–C). Where corals died, their tissue area became zero and they were not continued in the plots. Figure 2. Open in new tabDownload slide Mean (± s.e. confidence belts) percentage cover of live tissue on colonies of Acropora spp., Pocillopora verrucosa, and massive Porites spp. in the fixed portions of three flumes incubated at different pCO2 treatments (A–C). Where corals died, their tissue area became zero and they were not continued in the plots. When displayed using two-dimensional ordination, the sizes of the colonies in each treatment described multivariate community structures differing among months and the ordinations were similar among treatments (Figure 3). However, based on Spearman rank correlations comparing resemblance matrices for the communities created by Acropora spp. and P. verrucosa under the two elevated pCO2 treatments vs. the control (i.e. 396 µatm), community dynamics differed from that recorded under 396 µatm for both the 782 µatm (ρ = 0.606, pperm = 0.001) and 1434 µatm (ρ = 908, pperm = 0.001) treatments. These effects are expressed in the uniform rate at which monthly community analyses diverged over time at 396 µatm pCO2, the four clusters of monthly community structures at 782 µatm, and the tighter clustering of the first 10 months of community structure at 1434 µatm (Figure 3). Figure 3. Open in new tabDownload slide nMDS of coral reef communities composed of Acropora spp. and Pocillopora verrucosa in three flumes incubated under different pCO2 treatments. Analyses are based on the area of individual corals in the fixed portion of the community for each month, with missing values interpolated by linear regression, and areas set to zero if colonies died. Values were log(x + 1) transformed, and dissimilarities calculated as Bray–Curtis values. Circles are scaled to the summed cover of Acropora spp. and P. verrucosa, and are sequentially linked by a vector from the start to the end of the study. Figure 3. Open in new tabDownload slide nMDS of coral reef communities composed of Acropora spp. and Pocillopora verrucosa in three flumes incubated under different pCO2 treatments. Analyses are based on the area of individual corals in the fixed portion of the community for each month, with missing values interpolated by linear regression, and areas set to zero if colonies died. Values were log(x + 1) transformed, and dissimilarities calculated as Bray–Curtis values. Circles are scaled to the summed cover of Acropora spp. and P. verrucosa, and are sequentially linked by a vector from the start to the end of the study. The physical and chemical conditions created over the year within the treatments (Supplementary Table S1) were tested for their ability to explain variation in the biotic data using Spearman rank correlations after screening standardized abiotic data for co-linearity. These analyses were completed for December 2017–September 2018, which corresponded to the 10-month period with a complete set of physical and chemical data in the flumes; the analysis employed seven metrics as predictor variables (Supplementary Table S1). The change in community structure that occurred in the flume maintained at 1434 µatm pCO2 was best explained by MMT (pperm = 0.003), although significant explanatory power also was provided by MMT and the slope of changing mean daily temperature prior to sampling (SMDT) (pperm < 0.050), and MMT, SMDT and the integrated light intensity on the day of measurement (pperm < 0.050); no other combination of variables was significant (pperm > 0.050). Variation in community structure in the 396- and 782 µatm flume was not explained by any combination of abiotic drivers (p-perm ≥ 0.050) (Supplementary Table S1). Discussion The present study describes a year-long experiment in which fore reef communities with ecological relevance to co-located ecological time-series analyses (Edmunds 2019) were exposed to ambient and high pCO2 under seasonally varying environmental conditions. The outcome demonstrates the effects of high pCO2 on organismic net calcification (Gnet) that reflected largely the inhibitory effects anticipated from previous research (Chan and Connelly 2012; Comeau et al., 2016; Kornder et al., 2018). These effects, however, did not translate into modified rates at which the percentage cover of coral colonies changed over time. Instead, the percentage cover of corals in each flume was strongly affected by common temporal drivers of variation in coral growth, including temperature and light (Knutson et al., 1972; Lough and Barnes, 2000; Pratchett et al., 2015; Scheufen et al., 2017). While the outcomes of the incubations were similar among treatments in terms of final coral cover, the effects of high CO2 were evident in statistically significant, but small, shifts in the extent to which community structure differed between months in each flume (Figure 3). As we have described for a back reef community incubated for a year in the same flumes (Edmunds et al., 2019a), the present results for a community from a very different habitat and exposed to different environmental conditions are broadly similar in terms of the response of coral community structure to high pCO2. The present results for a fore reef community suggest that the effects of OA on coral community structure, and the contribution of these effects to the coral reef crisis (sensu Hoegh-Guldberg et al., 2007; Hughes et al., 2010), could remain subtle for years to come. While the generality of this conclusion is limited by the year duration of the experiment, longer duration experiments probably will remain intractable in remote tropical locations. To lengthen the effective exposure times in studies of the effects of OA on coral communities, future work should consider approaches differing from the one employed herein, for example utilizing greater use of predictive modelling (e.g. Evenhuis et al., 2014). The present study differs from most other studies addressing the effects of OA on corals and coral reefs because we focused on ecological effects in terms of a state variable (i.e. coral cover) that is at the forefront of quantifying the coral reef crisis (e.g. Bruno and Selig, 1970; De'ath et al., 2012; Jackson et al., 2014). Two decades of research on the effects of OA on the net calcification of corals, other coral reef organisms, and reef communities have revealed the implications of this stressor for the persistence of coral reefs as a calcified ecosystem (e.g. Silverman et al., 2009; Eyre et al., 2018). While the focus of this rapid advancement of knowledge is consistent with the importance of calcification for the functional significance of coral reefs (Anderson and Gledhill, 2013; Pratchett et al., 2015), and the role of scleractinians as a foundation taxon (Jones et al., 1994), it has left ecological processes overlooked as targets of OA effects. This has created a disconnect between ecological analyses of coral reefs through which the coral reef crisis has been recorded, largely through the reduction in coral cover (Bruno and Selig, 2007; De'ath et al., 2012; Jackson et al., 2014), and physiological and biogeochemical analyses through which the impairment of calcification by OA has been reported (Doney et al., 2009; Anderson and Gledhill, 2013). While it is likely that calcification of most corals and calcified algae will be depressed in a more acidic future, and that reefs will transition into net dissolution within the current century due to declining seawater pH (Silverman et al., 2009; Eyre et al., 2018), it is unclear how these trends will affect the occupation of planar space by corals. Within the research community investigating the effects of OA on marine systems, there is a growing recognition of the importance of addressing the effects on ecological processes such as primary production, the energetic costs of consumers, and species interactions (Gaylord et al., 2015). Studies of the effects of OA on ecological processes shaping coral communities is in its infancy, in part because many of these processes unfold over time scales associated with the generation time of the organism involved. The effects of OA on coral reproduction, for example will be expressed only over time scales embracing annual cycles of gametogenesis, spawning, and recruitment (Albright et al., 2010; Albright and Langdon, 2011; Fabricius et al., 2017). The effects on coral–coral competition will appear only if corals have the time to grow towards one another, interact, and develop a dominant and subordinate hierarchy (Horwitz et al., 2017). Despite these experimental difficulties, improved technical competence in conducting OA experiments and the discovery of “natural laboratories” where seawater pCO2 is normally high (Fabricius et al., 2011; Inoue et al., 2013) are starting to reveal how coral community ecology might be affected by OA. For instance, coral recruits show depressed skeletogenesis and deformities at high pCO2 (Foster et al., 2016), yet coral recruitment is unaffected by habitually high pCO2 on at least one reef adjacent to a CO2 seep (Oprandi et al., 2019). As recruits contribute to creating crowded coral communities, competition among adjacent colonies is affected by high pCO2, thereby modifying competitive hierarchies among the component corals (Horwitz et al., 2017). One potential outcome of these kinds of ecological effects is modified coral community structure. While the state variable most widely used to quantify coral communities, coral cover, remains understudied with respect to OA, there is evidence that it may be less strongly affected by OA than the net calcification of coral colonies. Support for this inference comes from two domains of investigation, one focusing on the ways in which mass deposition of CaCO3 translates into linear skeletal extension and the other focusing on measurements of the cover state variable itself. In terms of CaCO3 deposition and skeletal extension, there are several experimental studies showing that skeletal porosity increases as OA depresses mass deposition in corals, thereby allowing linear extension (and increase in coral colony cover) to be sustained (Fantazzini et al., 2015; Mollica et al., 2018), presumably at the expense of skeletal strength (Hennige et al., 2015). Studies of coral communities at some natural CO2 seeps show that the cover of scleractinians can be maintained (Fabricius et al., 2011), but not always (Inoue et al., 2013). Our previous work with back reef communities in Móorea shows that coral cover can persist while the carbonate framework experiences dissolution, at least over a year (Edmunds et al., 2019a). The present study provides a second example from a very different reef habitat of a similar outcome. A limitation of the present study is that many of the P. verrucosa and massive Porites in the flumes shrank in cover or died. As coral cover at 17-m depth on the fore reef of Mo’orea increased over the same period (Edmunds, 2019), it is reasonable to conclude that the declines in coral cover in the flumes were associated with the ex situ environment. Corals thrived in the same flumes when a similar experiment was completed with back reef communities over 2015–2016 (Comeau et al., 2015), and also when fore reef communities from 17-m depth were grown for 7 weeks from August to October 2014 (Comeau et al., 2016). Seawater temperature during the present study was slightly warmer in some months compared with the same months in the previous studies (i.e. 0.5°C warmer in January 2018 and 1.6°C warmer in September 2018 compared with the same months in 2016 (Edmunds et al., 2019a; see also Edmunds 2017), but none of the MMTs during 2017–2018 exceeded the bleaching threshold for Mo’orea (30°C, https://coralreefwatch.noaa.gov/). Since seawater temperatures during the present year-long experiment were not extreme relative to the local bleaching threshold, it is unlikely that thermal stress was the cause of coral mortality. While the reasons for the coral mortality in our flumes remain unknown, and despite pumping seawater from 14-m depth in Cook’s Bay, we cannot exclude the possibility that the coral mortality was associated with impaired seawater quality, or possibly, the negative effects of short exposure to high PFDs (e.g. Baker 2001) when the lids of the flumes were removed for cleaning. Previously, Comeau et al. (2016) used a 7-week experiment in the austral spring to show that 24-h net calcification of fore reef communities from 17-m depth in Moorea was reduced by 45% at 1176 vs. 401 µatm pCO2. This effect was caused by a 31% depression of daytime net community calcification and a 76% reduction in night-time net community calcification that was attributed to dissolution within reef rock (Comeau et al., 2016). The authors did not address how coral cover was affected by the treatments (because measurable planar growth was not expected over 7 weeks). However, with the similarity of physical and chemical conditions between the treatments in Comeau et al. (2016) and the present study, and the steps that were taken to make the fore reef communities similar in both studies (described in “Material and Methods” section), it is reasonable to expect that community calcification would be similarly depressed by OA in the present study, again, through dissolution within reef rock. Preliminary analyses (i.e. for November 2017–June 2018) for another aspect of the present study shows an 83% reduction in 24 h net calcification at 1434 vs. 396 µatm pCO2, with the greater effect (cf. Comeau et al., 2016) likely reflecting the decline in coral cover over time in the present study. Against this backdrop, it is striking that the effect of OA on community structure in the present study was confined to subtle shifts in the rate at which multivariate community structure composed of Acropora spp. and P. verrucosa separated among months over a year. Elucidating the causes of these effects was beyond the scope of the present study and, indeed, analyses of the association between multivariate community structure and multivariate environmental conditions proved equivocal in identifying dominant drivers of the changes in coral community structure. Given the seasonality in environmental conditions captured by the present experiment (Supplementary Table S1), and the growing understanding of interactive effects among OA and other environmental conditions in mediating coral performance (e.g. Chan et al., 2016; Langdon et al., 2018), it is likely that such effects were responsible for the shift in the tempo of variation in community structure in the present study. The present study contributes to a small number of studies that are experimentally addressing the effects of OA on coral communities over time scales that have ecological relevance for reef corals. Such experiments pose unique logistical challenges for the research community, and as the present study reveals, they can capture unexpected effects such as shrinkage and mortality of the organisms under investigation. Nevertheless, the outcome of such experiments can still be insightful when controls reveal the same trends and, as the present study shows (as well as our previous study, Edmunds et al., 2019a), can reveal findings unobtainable within the scale of short experiments (i.e. lasting weeks to months) conducting in small tanks. The present analysis reveals the extent to which the effects of OA on the planar coral cover (and coral community structure derived therefrom) can be undetected over a year. Further empirical or modelling approaches will be required to extend the time horizon towards the generation time of coral colonies (i.e. years–decades), or the time required for reefs to recover from disturbances in the Anthropocene (i.e. decades–centuries). Acknowledgements We thank the staff of the Richard B. Gump South Pacific Research Station for supporting our research. G. Srednick, S. Ginther, and S. Merolla provided key technical field support, and the year of flume operation was made possible with the help of A. Potter, A. Widrick, A. Isaac, J. Serrano, L. Perng, and B. Shayka. Research was completed under permits issued by the French Polynesian Government (Délégation à la Recherche) and the Haut-Commissariat de la République en Polynésie Française (DTRT) (Protocole d’Accueil 2017–2018). This is contribution number 299 of the Marine Biology Program of CSUN. Funding This research was funded by the US National Science Foundation (OCE 14-15268, 12-36905). References Abramoff M. D. , Magalhães P. J. , Ram S. J. 2004 . Image processing with ImageJ . Biophotonics International , 11 : 36 – 42 . OpenURL Placeholder Text WorldCat Albright R. , Takeshita Y. , Koweek D. A. , Ninokawa A. , Wolfe K. , Rivlin T. , Nebuchina Y. , et al. . 2018 . Carbon dioxide addition to coral reef waters suppresses net community calcification . Nature , 555 : 516 – 519 . Google Scholar Crossref Search ADS PubMed WorldCat Albright R. , Langdon C. 2011 . Ocean acidification impacts multiple early life history processes of the Caribbean coral Porites astreoides . Global Change Biology , 17 : 2478 – 2487 . Google Scholar Crossref Search ADS WorldCat Albright R. , Mason B. , Miller M. , Langdon C. 2010 . Ocean acidification compromises recruitment success of the threatened Caribbean coral Acropora palmata . Proceedings of the National Academy of Sciences of the United States of America , 107 : 20400 – 20404 . Google Scholar Crossref Search ADS PubMed WorldCat Anderson A. J. , Gledhill D. 2013 . Ocean acidification and coral reefs: effects on breakdown, dissolution, and net ecosystem calcification . Annual Review of Marine Science , 5 : 321 – 348 . Google Scholar Crossref Search ADS PubMed WorldCat Baker A. C. 2001 . Reef corals bleach to survive change . Nature , 411 : 765 – 766 . Google Scholar Crossref Search ADS PubMed WorldCat Bellwood D. R. , Hughes T. P. , Folke C. , Nyström M. 2004 . Confronting the coral reef crisis . Nature , 429 : 827 – 833 . Google Scholar Crossref Search ADS PubMed WorldCat Broecker W. S. 1975 . Climate change: are we on the brink of a pronounced global warming? Science , 189 : 460 – 463 . Google Scholar Crossref Search ADS PubMed WorldCat Bruno J. F. , Selig E. R. 2007 . Regional decline of coral cover in the Indo-Pacific: timing, extent, subregional comparisons . PLoS One , 8 : e711 . Google Scholar Crossref Search ADS WorldCat Carpenter R. 2018 . MCR LTER: coral reef: long-term population and community dynamics: benthic algae and other community components, ongoing since 2005. knb-lter-mcr.8.30, doi: 10.6073/pasta/9041bae6c2a2c2c61eb2d48684f7ac1e. Carpenter R. C. , Lantz C. A. , Shaw E. , Edmunds P. J. 2018 . Responses of coral reef community metabolism in flumes to ocean acidification . Marine Biology , 165 : 66 . Google Scholar Crossref Search ADS WorldCat Chan, N. C. S., Wangpraseurt, D., Kühl, M. and Connelly, S. R. 2016. Flow and coral morphology control coral surface pH: Implications for the effects of ocean acidification. Frontiers in Marine Science, 3: 10. Chan N. C. S. , Connolly S. R. 2012 . Sensitivity of coral calcification to ocean acidification: a meta-analysis . Global Change Biology , 19 : 282 – 90 . Google Scholar Crossref Search ADS PubMed WorldCat Clarke K. R. , Gorley R. N. 2006 . PRIMER V6: User Manual-Tutorial . PRIMER-E , Plymouth . 190 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Comeau S. , Carpenter R. C. , Lantz C. A. , Edmunds P. J. 2015 . Ocean acidification accelerates dissolution of experimental coral reef communities . Biogeosciences , 12 : 365 – 372 . Google Scholar Crossref Search ADS WorldCat Comeau S. , Lantz C. A. , Edmunds P. J. , Carpenter R. C. 2016 . Framework of barrier reefs threatened by ocean acidification . Global Change Biology , 22 : 1225 – 1234 . 10.1111/gcb.13023 . Google Scholar Crossref Search ADS PubMed WorldCat Crossref Connell J. H. 1973 . Population ecology of reef-building corals. In Biology and Geology of Coral Reefs Volume II: Biology 1 , pp. 205 – 244 . Ed. By Jones O. A. , Endean R. . Academic Press, New York and London . 464 pp. Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Cornwall C. E. , Hurd C. L. 2015 . Experimental design in ocean acidification research: problems and solutions . ICES Journal of Marine Science , 73 : 572 – 581 . Google Scholar Crossref Search ADS WorldCat De’ath G. , Fabricius K. E. , Sweatman H. , Puotinen M. 2012 . The 27-year decline of coral cover on the Great Barrier Reef and its causes . Proceedings of the National Academy of Sciences of the United States of America , 109 : 17995 – 17999 . Google Scholar Crossref Search ADS PubMed WorldCat Dickson A. G. , Wesolowski D. J. , Palmer D. A. , Mesmer R. E. 1990 . Dissociation constant of bisulfate ion in aqueous sodium chloride solution to 250°C . Journal of Physical Chemistry , 94 : 7978 – 7985 . Google Scholar Crossref Search ADS WorldCat Dickson A. G. , Sabine C. L. , Christian J. R. (Eds) 2007 . Guide to Best Practices for Ocean CO2 Measurements . PICES Special Publication 3, North Pacific Marine Science Organization , Sidney . 191 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Done T. J. 1983 . Coral zonation, its nature and significance. In Perspectives on Coral Reefs , pp. 107 – 147 . Ed. by Barnes D. J. , Clouston B. . Australian Institute of Marine Science , Australia . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Doney S. C. , Fabry V. J. , Feely R. A. , Kleypas J. A. 2009 . Ocean acidification: the other CO2 problem . Annual Review of Marine Science , 1 : 169 – 192 . Google Scholar Crossref Search ADS PubMed WorldCat Donner S. D. , Heron S. F. , Skirving W. J. 2018 . Future scenarios: a review of modeling efforts to predict the future of coral reefs in an era of climate change. In Coral Bleaching Patterns, Processes, Causes and Consequences , 2nd edn, pp. 325 – 341 . Ed. by van Oppen M. J. H. , Lough J. M. . Ecological Studies (Analysis and Synthesis), Springer , Cham . 356 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Doo S. S. , Edmunds P. J. , Carpenter R. C. 2019 . Ocean acidification effects on in situ coral reef metabolism . Scientific Reports , 9 : 12067 . Google Scholar Crossref Search ADS PubMed WorldCat Dove S. G. , Kline D. I. , Pantos O. , Angly F. E. , Tyson G. W. , Hoegh-Guldberg O. 2013 . Future reef decalcification under a business-as-usual CO2 emission scenario . Proceedings of the National Academy of Sciences of the United States of America , 110 : 15342 – 15347 . Google Scholar Crossref Search ADS PubMed WorldCat Edmunds P. J. 2015 . A quarter-century demographic analysis of the Caribbean coral, Orbicella annularis, and projections of population size over the next century . Limnology and Oceanography , 60 : 840 – 855 . Google Scholar Crossref Search ADS WorldCat Edmunds P. J. 2017 . Unusually high coral recruitment during the 2016 El Niño in Mo’orea, French Polynesia . PLoS One , 12 : e0185167 . Google Scholar Crossref Search ADS PubMed WorldCat Edmunds P. J. 2018 . Implications of high rates of sexual recruitment in driving rapid reef recovery in Mo’orea, French Polynesia . Scientific Reports , 8 : 16615 . Google Scholar Crossref Search ADS PubMed WorldCat Edmunds P. 2019 . MCR LTER: coral reef: long-term population and community dynamics: corals, ongoing since 2005. knb-lter-mcr.4.37, doi:10.6073/pasta/721fdcde85b630d1b4548b677d9b1d86. Edmunds P. J. , Adam T. C. , Baker A. C. , Doo S. S. , Glynn P. W. , Manzello D. P. , Silbiger N. J. et al. . 2019 . Why more comparative approaches are required in time-series analyses of coral reef ecosystems . Marine Ecology Progress Series , 608 : 297 – 306 . Google Scholar Crossref Search ADS WorldCat Edmunds P. J. , Burgess S. C. 2016 . Size-dependent physiological responses of the branching coral Pocillopora verrucosa to elevated temperature and pCO2 . Journal of Experimental Biology , 219 : 3896 – 3906 . Google Scholar Crossref Search ADS PubMed WorldCat Edmunds P. J. , Doo S. S. , Carpenter R. C. 2019 . Changes in coral reef community structure in response to year-long incubations under contrasting pCO2 regimes . Marine Biology , 166 : 94 – 105 . Google Scholar Crossref Search ADS WorldCat Edmunds P. J. , Leichter J. J. , Johnston E. C. , Tong E. J. , Toonen R. J. 2016 . Ecological and genetic variation in reef-building corals on four Society Islands . Limnology and Oceanography , 61 : 543 – 557 . Google Scholar Crossref Search ADS WorldCat Evensen N. R. , Edmunds P. J. 2017 . Conspecific aggregations mitigate the effects of ocean acidification on calcification of the coral Pocillopora verrucosa . Journal of Experimental Biology , 220 : 1097 – 1105 . Google Scholar Crossref Search ADS PubMed WorldCat Evenhuis C. , Lenton A. , Cantin N. E. , Lough J. N. 2014 . Modeling coral calcification accounting for the impacts of coral bleaching and ocean acidification . Biogeosciences , 11 : 187 – 249 . Google Scholar Crossref Search ADS WorldCat Eyre B. D. , Cyronak T. , Drupp P. , Heinen De Carlo E. , Sachs J. P. , Andersson A. J. 2018 . Coral reefs will transition to net dissolving before end of century . Science , 359 : 908 – 911 . Google Scholar Crossref Search ADS PubMed WorldCat Fabricius K. E. , Langdon C. , Uthicke S. , Humphrey C. , Noonan S. , De’ath G. , Okazaki R. et al. . 2011 . Losers and winners in coral reefs acclimatized to elevated carbon dioxide concentrations . Nature Climate Change , 1 : 165 – 169 . Google Scholar Crossref Search ADS WorldCat Fabricius K. E. , Noonan S. H. C. , Abrego D. , Harrington L. , De’ath G. 2017 . Low recruitment due to altered settlement substrata as primary constraint for coral communities under ocean acidification. Proceedings of the Royal Society B Biological Sciences. 284: doi: 10.1098/rspb.2017.1536. Fantazzini P. , Mengoli S. , Pasquini L. , Bortolotti V. , Brizi L. , Mariani M. , Di Giosia M. et al. . 2015 . Gains and losses of coral skeletal porosity changes with ocean acidification acclimation . Nature Communications , 6 : doi: 10.1038/ncomms8785. OpenURL Placeholder Text WorldCat Foster T. , Falter J. L. , McCulloch M. T. , Clode P. L. 2016 . Ocean acidification causes structural deformities in juvenile coral skeletons . Science Advances , 2 : e1501130 . Google Scholar Crossref Search ADS PubMed WorldCat Gaylord B. , Kroeker K. J. , Sunday J. M. , Anderson K. M. , Barry J. P. , Brown N. E. , Connell S. D. et al. . 2015 . Ocean acidification through the lens of ecological theory . Ecology , 96 : 3 – 15 . Google Scholar Crossref Search ADS PubMed WorldCat Hench J. L. , Leichter J. J. , Monismith S. G. 2008 . Episodic circulation and exchange in a wave-driven coral reef and lagoon system . Limnology and Oceanography , 53 : 2681 – 1694 . Google Scholar Crossref Search ADS WorldCat Hennige S. J. , Wicks L. C. , Kamenos N. A. , Perna G. , Findlay H. S. , Roberts J. M. 2015 . Hidden impacts of ocean acidification to live and dead coral framework. Proceedings of the Royal Society B Biological Sciences, 282: doi: 10.1098/rspd.2015.0990. Hoegh-Guldberg O. , Mumby P. J. , Hooten A. J. , Steneck R. S. , Greenfield P. , Gomez E. , Harvell C. D. et al. . 2007 . Coral reefs under rapid climate change and ocean acidification . Science , 318 : 1737 – 1742 . Google Scholar Crossref Search ADS PubMed WorldCat Hofmann G. E. , Smith J. E. , Johnson K. S. , Send E. , Levin L. A. , Micheli F. , Paytan A. et al. . 2011 . High frequency dynamics of ocean pH: a multi-ecosystem comparison . PLoS One , 6 : e28983 . Google Scholar Crossref Search ADS PubMed WorldCat Horwitz R. , Hoogenboom M. O. , Fine M. 2017 . Spatial competition dynamics between reef corals under ocean acidification . Scientific Reports , 7 : doi: 10.1038/srep40288. OpenURL Placeholder Text WorldCat Hurlbert S. H. 1984 . Pseudoreplication and the design of ecological field experiments . Ecological Monographs , 54 : 187 – 211 . Google Scholar Crossref Search ADS WorldCat Hughes T., P. , Graham N. A. J. , Jackson J. B. C. , Mumby P. J. , Steneck R. S. 2010 . Rising to the challenge of sustaining coral reef resilience . Trends in Ecology & Evolution , 25 : 633 – 642 . Google Scholar Crossref Search ADS PubMed WorldCat Inoue S. , Kayanne H. , Yamamoto S. , Kurihara H. 2013 . Spatial community shift from hard to soft corals in acidified water . Nature Climate Change , 3 : 683 – 687 . Google Scholar Crossref Search ADS WorldCat IPCC. 2014. Climate Change 2014: Synthesis Report. In Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Ed. by Core Writing Team, R. K., Pachauri, L. A. Meyer. IPCC, Geneva, pp. 151 . Jackson J. B. C. , Donovan M. K. , Cramer K. L. , Lam V. V. (Eds) 2014 . Status and trends of Caribbean coral reefs: 1970-2012. In Global Coral Reef Monitoring Network. IUCN, Gland. 304 pp. Jokiel P. L. , Rodgers K. S. , Kuffner I. B. , Andersson A. J. , Cox E. F. , Mackenzie F. T. 2008 . Ocean acidification and calcifying reef organisms: a MESOCOSM investigation . Coral Reefs , 27 : 473 – 483 . Google Scholar Crossref Search ADS WorldCat Jones C. G. , Lawton J. H. , Shachak M. 1994 . Organisms as ecosystem engineers . Oikos , 69 : 373 – 386 . Google Scholar Crossref Search ADS WorldCat Kleypas J. A. , Buddemeier R. W. , Archer D. , Gattuso J. , Langdon C. , Opdyke B. N. 1999 . Geochemical consequences of increased atmospheric carbon dioxide on coral reefs . Science , 284 : 118 – 120 . Google Scholar Crossref Search ADS PubMed WorldCat Kline D. I. , Teneva L. , Schneider K. , Miard T. , Chai A. , Marker M. , Headley K. et al. . 2012 . A short-term in situ CO2 enrichment experiment on Heron Island (GBR) . Scientific Reports , 2 : doi: 10.1038/srep00413. OpenURL Placeholder Text WorldCat Knutson D. W. , Buddemeier R. W. , Smith S. V. 1972 . Coral chronometers: seasonal growth bands in reef corals . Science , 177 : 270 – 272 . Google Scholar Crossref Search ADS PubMed WorldCat Kornder N. A. , Riegl B. M. , Figueiredo J. 2018 . Thresholds and drivers of coral calcification responses to climate change . Global Change Biology , 24 : 5084 – 5095 . Google Scholar Crossref Search ADS PubMed WorldCat Kroeker K. J. , Kordas R. L. , Crim R. , Hendriks I. E. , Ramajo L. , Singh G. S. , Duarte C. M. , et al. . 2013 . Impacts of ocean acidification on marine organisms: quantifying sensitivities and interaction with warming . Global Change Biology , 19 : 1884 – 1896 . Google Scholar Crossref Search ADS PubMed WorldCat Kroeker K. J. , Kordas R. L. , Crim R. N. , Singh G. G. 2010 . Meta-analysis reveals negative yet variable effects of ocean acidification on marine organisms . Ecology Letters , 13 : 1419 – 1434 . Google Scholar Crossref Search ADS PubMed WorldCat Langdon C. , Albright R. , Baker A. C. , Jones P. 2018 . Two threatened Caribbean coral species have contrasting responses to combined temperature and acidification stress . Limnology and Oceanography , 63 : 2450 – 2464 . Google Scholar Crossref Search ADS WorldCat Langdon C. , Takahashi T. , Sweeney C. , Chipman D. , Goddard J. , Marubini F. , Aceves H. et al. . 2000 . Effect of calcium carbonate saturation state on the calcification rate of an experimental coral reef . Global Biogeochemical Cycles , 14 : 639 – 654 . Google Scholar Crossref Search ADS WorldCat Lavigne H. and Gattuso J. P. 2013. Seacarb, seawater carbonate chemistry with R. R package version 2.4.10. Available at: https://CRAN.R-project.org/package=seacarb (last accessed 11 July 2014). Lough J. M. , Barnes D. J. 2000 . Environmental controls on growth of the massive coral Porites . Journal of Experimental Marine Biology and Ecology , 245 : 225 – 243 . Google Scholar Crossref Search ADS PubMed WorldCat Loya Y. 1972 . Community structure and species diversity of hermatypic corals at Eilat, Red Sea . Marine Biology , 13 : 100 – 123 . Google Scholar Crossref Search ADS WorldCat Lueker T. J. , Dickson A. G. , Keeling C. D. 2000 . Ocean pCO2 calculated from dissolved inorganic carbon, alkalinity, and equations for K-1 and K-2: validation based on laboratory measurements of CO2 in gas and seawater at equilibrium . Marine Chemistry , 70 : 105 – 119 . Google Scholar Crossref Search ADS WorldCat Mollica N. R. , Guo W. , Cohen A. L. , Huang K. , Foster G. L. , Donald H. K. , Solow A. R. 2018 . Ocean acidification affects coral growth by reducing skeletal density . Proceedings of the National Academy of Sciences of the United States of America , 115 : 1754 – 1759 . Google Scholar Crossref Search ADS PubMed WorldCat Oprandi A. , Montefalcone M. , Morri C. , Benelli F. , Bianchi C. N. 2019 . Water circulation, and not ocean acidification, affects coral recruitment and survival at shallow hydrothermal vents . Estuarine, Coastal and Shelf Science , 217 : 158 – 164 . Google Scholar Crossref Search ADS WorldCat Perez F. F. , Fraga F. 1987 . The pH measurements in seawater on the NBS scale . Marine Chemistry , 21 : 315 – 327 . Google Scholar Crossref Search ADS WorldCat Pratchett M. S. , Anderson K. D. , Hoogenboom M. O. , Widman E. , Baird A. H. , Pandol J. M. , Edmunds P. J. et al. . 2015 . Spatial, temporal and taxonomic variation in coral growth-implications for the structure and function of coral reef ecosystems. In Oceanography and Marine Biology an Annual Review , Volume 53 , 1st edn, pp. 215 – 295 . Ed. by Hughes R. N. , Hughes D. J. , Smith I. P. , Dale A. C. . CRC Press , Boca Raton . 366 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rippe J. P. , Baumann J. H. , De Leener D. N. , Aichelman H. E. , Friedlander E. B. , Davies S. W. , Castillo K. D. 2018 . Corals sustain growth but not skeletal density across the Florida Keys Reef Tract despite ongoing warming . Global Change Biology , 24 : 5205 – 5217 . Google Scholar Crossref Search ADS PubMed WorldCat Scheufen T. , Krämer W. E. , Iglesias-Prieto R. , Enríquez S. 2017 . Seasonal variation modulates coral sensibility to heat-stress and explains annual changes in coral productivity . Scientific Reports , 7 : doi: 10.1038/s41598-017-04927-8. OpenURL Placeholder Text WorldCat Silverman J. , Lazar B. , Cao L. , Caldeira K. , Erez J. 2009 . Coral reefs may start dissolving when atmospheric CO2 doubles . Geophysical Research Letters , 36 : doi: 10.1029/2008GL036282. OpenURL Placeholder Text WorldCat Spencer Davies P. S. 1989 . Short-term growth measurements of corals using an accurate buoyant weighing technique . Marine Biology , 101 : 389 – 395 . Google Scholar Crossref Search ADS WorldCat Vecsei A. 2001 . Fore-reef carbonate production: development of a regional census-based method and first estimate . Palaeogeography, Palaeoclimatology, Palaeoecology , 175 : 185 – 200 . Google Scholar Crossref Search ADS WorldCat Vecsei A. 2004 . A new estimate of global reefal carbonate production including the fore reef . Global Planetary Change , 43 : 1 – 18 . Google Scholar Crossref Search ADS WorldCat Veron J. E. N. 2000 . Corals of the World . Ed. by Stafford-Smith M. . Australian Institute of Marine Science , Townsville . 489 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC van Hooidonk R. , Maynard J. A. , Manzello D. , Planes S. 2014 . Opposite latitudinal gradients in projected ocean acidification and bleaching impacts on coral reefs . Global Change Biology , 20 : 103 – 112 . Google Scholar Crossref Search ADS PubMed WorldCat Washburn L. 2018 . MCR LTER: Coral Reef: Ocean Currents and Biogeochemistry: salinity, temperature and current at CTD and ADCP mooring FOR01 from 2004 ongoing. knb-lter-mcr.30.32 doi: 10.6073/pasta/2f4279e07911b128f11aaf7503eb5c02. © International Council for the Exploration of the Sea 2020. All rights reserved. For permissions, please email: [email protected] This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Morse, Molly, R;Kerr, Lisa, A;Galuardi,, Benjamin;Cadrin, Steven, X

doi: 10.1093/icesjms/fsaa082pmid: N/A

Abstract Accounting for movement and mixing in stock assessment is important for managing sustainable fisheries, particularly for highly migratory species. However, many fisheries management approaches continue to use single-stock, single-area models to assess mixed-population stocks that are known to have complex movement dynamics. We evaluated a single-stock, single-area stock assessment model’s performance on fishery pseudodata generated using a spatially complex operating model that incorporates movement and mixing of simulated Atlantic bluefin tuna-like populations. Structural model misspecification produced positively biased perceptions of size and productivity of the smaller western population, based on supplement by the larger eastern population, and negatively biased perceptions of the size and productivity of the eastern population due to net movement of fish out of the eastern stock area. This bias could lead to unintended overexploitation of the smaller western population and potential for foregone yield of the larger eastern population. Our findings provide a greater understanding of the effects of movement and mixing on single-stock, single-area model-based management approaches and emphasize the importance of explicitly considering these dynamics in ensuring the sustainability of highly migratory species like Atlantic bluefin tuna. Introduction Complexities in the spatial distribution and movement patterns of exploited fish populations pose challenges to stock assessment modelling and fishery management (Cadrin and Secor, 2009). The implications of these complexities are magnified for valuable, highly migratory species targeted by several nations with diverse interests, priorities, and economies (Webster, 2009; Pons et al., 2018). Understanding the patterns of movement and mixing for such highly migratory species and accounting for them in stock assessment and fishery management are important to ensuring the sustainability of these fisheries. Spatial approaches to stock assessment can improve estimation performance compared to those that ignore spatial structure (Punt, 2019b). Simulation testing is a particularly powerful tool for evaluating stock assessment model performance, and several studies have used it to test the consequences of spatial structure and connectivity among fish populations on stock assessment outcomes (Guan et al., 2013; Deroba et al., 2015; Szuwalski and Punt, 2015). Kerr and Goethel (2014) and Goethel et al. (2016) developed a general approach for constructing spatially explicit operating models to test and validate stock assessment models. In this approach, a spatial assessment model is fit to available data to obtain parameter estimates that are used to condition a spatially explicit operating model. The operating model can simulate various population scenarios, and the performance of stock assessment models and management procedures can be evaluated. Goethel et al. (2015) used this approach to demonstrate that a spatially explicit, tag-integrated model with reliable tagging data performed better than a closed-population model for a metapopulation of three yellowtail flounder-like stocks off New England. The Atlantic bluefin tuna (Thunnus thynnus) is a highly migratory species that inhabits the North Atlantic Ocean and Mediterranean Sea and exhibits complex intermixing of genetically distinct western and eastern populations (National Research Council (NRC), 1994; Fromentin and Powers, 2005). Otolith chemistry, tagging, and genetic data suggest that the western and eastern populations exhibit spawning site fidelity to the Gulf of Mexico and Mediterranean Sea, respectively (Nemerson et al., 2000; Block et al., 2005; Boustany et al., 2008; Rooker et al., 2014). These populations mix extensively in North Atlantic feeding grounds, but mixing varies across space, time, and demographic groups (Rooker et al., 2008a, b, 2014; Siskey et al., 2016; Morse et al., 2018; Puncher et al., 2018). However, stock assessment methods used to inform management of these fisheries do not account for this important aspect of stock structure. The International Commission for the Conservation of Atlantic Tunas (ICCAT) manages the Atlantic bluefin tuna fishery with a two-stock approach, which delineates the western and eastern stocks by a boundary at 45°W in the North Atlantic (ICCAT, 1981; Figure 1). ICCAT conducts separate western and eastern stock assessments assuming closed populations without mixing (i.e. “unit stocks”; Secor, 2014). In the past, stock assessment models developed to account for population mixing of Atlantic bluefin tuna (Butterworth and Punt, 1994; Porch, 1995, 2003; Porch et al., 1995; Taylor et al., 2011) did not consistently perform better than simpler models (Porch et al., 1998, 2001) and data on population mixing were limited (Porch, 2005). ICCAT has used this as justification to continue managing these fisheries based on separate western and eastern stock assessments. However, ICCAT has acknowledged the great uncertainty that mixing poses to stock assessment results (ICCAT Standing Committee on Research and Statistics (SCRS), 2017; ICCAT, 2018a). Figure 1. Open in new tabDownload slide Atlantic bluefin tuna stock management areas (Fromentin, 2006) and operating model geographic zones: Gulf of Mexico (zone 1), Gulf of St. Lawrence (zone 2), western Atlantic (zone 3), central Atlantic (zone 4), eastern Atlantic (zone 5), northeast Atlantic (zone 6), and Mediterranean Sea (zone 7; Kerr et al., 2017). Figure 1. Open in new tabDownload slide Atlantic bluefin tuna stock management areas (Fromentin, 2006) and operating model geographic zones: Gulf of Mexico (zone 1), Gulf of St. Lawrence (zone 2), western Atlantic (zone 3), central Atlantic (zone 4), eastern Atlantic (zone 5), northeast Atlantic (zone 6), and Mediterranean Sea (zone 7; Kerr et al., 2017). The recent availability of data sets to estimate movement and mixing rates (e.g. Boustany et al., 2008; Rooker et al., 2014; Siskey et al., 2016; Galuardi et al., 2014, 2018) makes Atlantic bluefin tuna an informative case study for testing the effects of mixing on stock assessment model performance. We employed a simulation testing framework to evaluate the performance of the single-stock, single-area virtual population analysis (VPA) stock assessment models developed by ICCAT on simulated mixed Atlantic bluefin tuna populations. We evaluated whether stock mixing would bias the results of these stock assessment models and impact fishery managers’ ability to make effective management decisions. The results from this study will inform future assessment and management of Atlantic bluefin tuna as well as contribute to the understanding of spatially complex fisheries and the implications of incorrect assumptions about spatial structure and population dynamics. Methods Operating model The analysis was based on the operating models developed by Kerr et al. (2017, 2018), which simulate two spawning populations: a western population originating in the Gulf of Mexico and an eastern population originating in the Mediterranean Sea. The operating model was deterministic and structured by age, fleet, season, and area, simulating movement of fish across seven geographic zones within two stock management areas (Figure 1). Seasonal fish movements among the seven zones were based on transfer rates derived from fishery-independent satellite tagging data (Galuardi et al., 2014). Movement was estimated independently for western juveniles (ages 1–8), western adults (ages 9–29), and eastern fish (all ages 1–29) to reflect current understanding of population and ontogenetic patterns in movement rates. Fish from one stock area could move to the other area, but only spawn in their natal area. Movement estimates were not constrained to require that fish return to their natal area during the spawning season, so reproductive contributions in the operating model were limited to mature fish that returned to the spawning grounds during the spawning season. We used the mean estimated movement rates as the base case movement scenario informing the operating model (Figure 2). To evaluate the effect of movement rate magnitude on estimation performance, we simulated an alternative movement scenario using the lower 25th percentile of movement rates from the range of estimates (Supplementary Figure S1). Figure 2. Open in new tabDownload slide Seasonal movement probability matrices for the base case movement scenario for (a) western juvenile fish (ages 1–8), (b) western adult fish (ages 9–29), and (c) eastern fish (all ages 1–29). Matrix values represent the proportion of fish moving from the starting geographic zones (y-axis) to the destination zones (x-axis; see Figure 1 for zones). Figure 2. Open in new tabDownload slide Seasonal movement probability matrices for the base case movement scenario for (a) western juvenile fish (ages 1–8), (b) western adult fish (ages 9–29), and (c) eastern fish (all ages 1–29). Matrix values represent the proportion of fish moving from the starting geographic zones (y-axis) to the destination zones (x-axis; see Figure 1 for zones). We revised the operating model from previous applications (Kerr et al., 2017, 2018) by conditioning it on estimates of recruitment, fishing mortality, natural mortality, spawning fraction, and growth from the 2017 ICCAT stock assessments of western and eastern Atlantic bluefin tuna (ICCAT, 2018a; Table 1 and Supplementary data A) [The spawning fraction is “the proportion of fish contributing to the spawning output of the population as a function of age”, which allows for the possibility that older mature fish spawn more frequently than younger mature fish (Porch and Hanke, 2018)]. The operating model was initialized using these estimates of abundance-at-age for the first modelled year (1974), and recruitment (age 1) and fishing mortality rates-at-age for all modelled years (1974–2015). Fishing mortality rates by age, year, geographic zone, and seasonal quarter were derived from stock-specific fishing mortality estimates, assumed quarterly proportions based on fishery patterns, and fleet-specific relative age-composition data (ICCAT, 2018a; eq. 4; Tables 2 and 3). Population and fishery dynamics were simulated over 42 years from 1974 to 2015 (eq. 5, 6; Tables 2 and 3). Table 1. Operating model specifications for western and eastern stocks and populations. Parameter . West . East . Stocks . Geographic zones in stock area (Figure 1) 1–3 4–7 Number of fleets 17 10 Populations Age classes Operating model Observation model 1–29 1–16+ Operating model Observation model 1-29 1-10+ Length-at-age Richards model (Table 2eq. 1) von Bertalanffy model (Table 2eq. 2) L1 L2 b A1 A2 K 33.0 270.6 −0.12 0 34 0.22 K L∞ a0 0.093 319 −0.97 Length-weight (Table 2 eq. 3) α β 1.77054E−5 3.00125 α β 3.50801E−5 2.87845 Natural mortality rate (annual)  Age 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16+  West 0.38 0.30 0.24 0.20 0.18 0.16 0.14 0.13 0.12 0.12 0.11 0.11 0.11 0.10 0.10 0.10  East 0.38 0.30 0.24 0.20 0.18 0.16 0.14 0.13 0.12 0.10 0.10 0.10 0.10 0.10 0.10 0.10 Spawning fraction  West 0 0 0 0 0 0.001 0.007 0.039 0.186 0.563 0.879 0.976 0.996 0.999 1 1  East 0 0 0.25 0.5 1 1 1 1 1 1 1 1 1 1 1 1 Parameter . West . East . Stocks . Geographic zones in stock area (Figure 1) 1–3 4–7 Number of fleets 17 10 Populations Age classes Operating model Observation model 1–29 1–16+ Operating model Observation model 1-29 1-10+ Length-at-age Richards model (Table 2eq. 1) von Bertalanffy model (Table 2eq. 2) L1 L2 b A1 A2 K 33.0 270.6 −0.12 0 34 0.22 K L∞ a0 0.093 319 −0.97 Length-weight (Table 2 eq. 3) α β 1.77054E−5 3.00125 α β 3.50801E−5 2.87845 Natural mortality rate (annual)  Age 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16+  West 0.38 0.30 0.24 0.20 0.18 0.16 0.14 0.13 0.12 0.12 0.11 0.11 0.11 0.10 0.10 0.10  East 0.38 0.30 0.24 0.20 0.18 0.16 0.14 0.13 0.12 0.10 0.10 0.10 0.10 0.10 0.10 0.10 Spawning fraction  West 0 0 0 0 0 0.001 0.007 0.039 0.186 0.563 0.879 0.976 0.996 0.999 1 1  East 0 0 0.25 0.5 1 1 1 1 1 1 1 1 1 1 1 1 All biological parameters are based on sources cited in ICCAT (2018a). Additional model specifications are documented in Supplementary data A. Open in new tab Table 1. Operating model specifications for western and eastern stocks and populations. Parameter . West . East . Stocks . Geographic zones in stock area (Figure 1) 1–3 4–7 Number of fleets 17 10 Populations Age classes Operating model Observation model 1–29 1–16+ Operating model Observation model 1-29 1-10+ Length-at-age Richards model (Table 2eq. 1) von Bertalanffy model (Table 2eq. 2) L1 L2 b A1 A2 K 33.0 270.6 −0.12 0 34 0.22 K L∞ a0 0.093 319 −0.97 Length-weight (Table 2 eq. 3) α β 1.77054E−5 3.00125 α β 3.50801E−5 2.87845 Natural mortality rate (annual)  Age 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16+  West 0.38 0.30 0.24 0.20 0.18 0.16 0.14 0.13 0.12 0.12 0.11 0.11 0.11 0.10 0.10 0.10  East 0.38 0.30 0.24 0.20 0.18 0.16 0.14 0.13 0.12 0.10 0.10 0.10 0.10 0.10 0.10 0.10 Spawning fraction  West 0 0 0 0 0 0.001 0.007 0.039 0.186 0.563 0.879 0.976 0.996 0.999 1 1  East 0 0 0.25 0.5 1 1 1 1 1 1 1 1 1 1 1 1 Parameter . West . East . Stocks . Geographic zones in stock area (Figure 1) 1–3 4–7 Number of fleets 17 10 Populations Age classes Operating model Observation model 1–29 1–16+ Operating model Observation model 1-29 1-10+ Length-at-age Richards model (Table 2eq. 1) von Bertalanffy model (Table 2eq. 2) L1 L2 b A1 A2 K 33.0 270.6 −0.12 0 34 0.22 K L∞ a0 0.093 319 −0.97 Length-weight (Table 2 eq. 3) α β 1.77054E−5 3.00125 α β 3.50801E−5 2.87845 Natural mortality rate (annual)  Age 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16+  West 0.38 0.30 0.24 0.20 0.18 0.16 0.14 0.13 0.12 0.12 0.11 0.11 0.11 0.10 0.10 0.10  East 0.38 0.30 0.24 0.20 0.18 0.16 0.14 0.13 0.12 0.10 0.10 0.10 0.10 0.10 0.10 0.10 Spawning fraction  West 0 0 0 0 0 0.001 0.007 0.039 0.186 0.563 0.879 0.976 0.996 0.999 1 1  East 0 0 0.25 0.5 1 1 1 1 1 1 1 1 1 1 1 1 All biological parameters are based on sources cited in ICCAT (2018a). Additional model specifications are documented in Supplementary data A. Open in new tab Table 2. Equations used in this study. Equation number . Operating model . 1 La,p=2=L1b+L2b-L1b1-e-Ka-A11-e-KA2-A11b 2 La,p=1=L∞1-e-Ka-a0 3 Wa=αLaβ 4 Fy,a,q,z=fy,a,sDz,qxy,a,z,g∑West z=1:3East z=4:7xy,a,z,g 5 Ny,a,q,z,p=Ny-1,a-1,q=4,z,pTa,q,z→z′,pe-(Fy,a,q,z+Ma-1,q,p) when q=1Ny,a,q-1,z,pTa,q,z→z′,pe-(Fy,a,q,z+Ma,q,p) when q=2:4 6 SSBy,q,z,p=∑a=129Ny,a,q,z,pBa,pWa,p1000 7 Ny,a,q,s=∑West z=1:3East z=4:7Ny,a,q,z,p 8 Ny,a,q,z=∑West p=2East p=1Ny,a,q,z,p 9 Cy,a,s=(∑West z=1:3East z=4:7∑q=14∑p=12Ny,a,q,z,p Fy,a,q,zFy,a,q,z+Ma,q,p1-e-Fy,a,q,z+Ma,q,p)eεy,a,s 10 Iy,g,s=∑West z=1:3East z=4:7∑a=129∑p=12Sa,gNy,a,q,z,pWa,pQgeεy,g,s 11 Ey,g=∑aAxy,a,z,giy,g,s 12 Xy,a,g,s=(∑West z=1:3East z=4:7∑q=14∑p=12Ny,a,q,z,pEy,gQgSa,gEy,gQgSa,g+Ma,q,p1-e-Ey,gQgSa,g+Ma,q,p)eεy,a,s Estimation model 13 Ny,a=Fy,a+My,aFy,a1-e-Fy,a+My,aCy,a 14 Fy-1,a-1=Fy-1,a-1+My-1,a-1Cy-1,a-1Ny,aeFy-1,a-1+My-1,a-1-1 15 Ny,A=Fy,A+My,AFy,A1-e-Fy,A+My,ACy,A 16 Ny-1,A=Fy-1,A+My-1,AFy-1,A1-e-Fy-1,A+My-1,ACy-1,A 17 Fy-1,A-1=Fy-1,A-1+My-1,A-1Cy-1,A-1eFy-1,A-1+My-1,A-1-1Ny,A-Ny-1,AeFy-1,A+My-1,A 18 I^y,g,s=Qy,g,s∑aAΔy,a,g,sSy,a,g,sWy,a,g,sNy,a,s Performance metrics 19 relative bias= medianymediannθestn,y-θtrueyθtruey×100% 20 absolute bias=medianymediannθestn,y-θtruey Equation number . Operating model . 1 La,p=2=L1b+L2b-L1b1-e-Ka-A11-e-KA2-A11b 2 La,p=1=L∞1-e-Ka-a0 3 Wa=αLaβ 4 Fy,a,q,z=fy,a,sDz,qxy,a,z,g∑West z=1:3East z=4:7xy,a,z,g 5 Ny,a,q,z,p=Ny-1,a-1,q=4,z,pTa,q,z→z′,pe-(Fy,a,q,z+Ma-1,q,p) when q=1Ny,a,q-1,z,pTa,q,z→z′,pe-(Fy,a,q,z+Ma,q,p) when q=2:4 6 SSBy,q,z,p=∑a=129Ny,a,q,z,pBa,pWa,p1000 7 Ny,a,q,s=∑West z=1:3East z=4:7Ny,a,q,z,p 8 Ny,a,q,z=∑West p=2East p=1Ny,a,q,z,p 9 Cy,a,s=(∑West z=1:3East z=4:7∑q=14∑p=12Ny,a,q,z,p Fy,a,q,zFy,a,q,z+Ma,q,p1-e-Fy,a,q,z+Ma,q,p)eεy,a,s 10 Iy,g,s=∑West z=1:3East z=4:7∑a=129∑p=12Sa,gNy,a,q,z,pWa,pQgeεy,g,s 11 Ey,g=∑aAxy,a,z,giy,g,s 12 Xy,a,g,s=(∑West z=1:3East z=4:7∑q=14∑p=12Ny,a,q,z,pEy,gQgSa,gEy,gQgSa,g+Ma,q,p1-e-Ey,gQgSa,g+Ma,q,p)eεy,a,s Estimation model 13 Ny,a=Fy,a+My,aFy,a1-e-Fy,a+My,aCy,a 14 Fy-1,a-1=Fy-1,a-1+My-1,a-1Cy-1,a-1Ny,aeFy-1,a-1+My-1,a-1-1 15 Ny,A=Fy,A+My,AFy,A1-e-Fy,A+My,ACy,A 16 Ny-1,A=Fy-1,A+My-1,AFy-1,A1-e-Fy-1,A+My-1,ACy-1,A 17 Fy-1,A-1=Fy-1,A-1+My-1,A-1Cy-1,A-1eFy-1,A-1+My-1,A-1-1Ny,A-Ny-1,AeFy-1,A+My-1,A 18 I^y,g,s=Qy,g,s∑aAΔy,a,g,sSy,a,g,sWy,a,g,sNy,a,s Performance metrics 19 relative bias= medianymediannθestn,y-θtrueyθtruey×100% 20 absolute bias=medianymediannθestn,y-θtruey Equation elements are defined in Tables 1 and 3. Estimation model equations are adapted from Porch (2003). Open in new tab Table 2. Equations used in this study. Equation number . Operating model . 1 La,p=2=L1b+L2b-L1b1-e-Ka-A11-e-KA2-A11b 2 La,p=1=L∞1-e-Ka-a0 3 Wa=αLaβ 4 Fy,a,q,z=fy,a,sDz,qxy,a,z,g∑West z=1:3East z=4:7xy,a,z,g 5 Ny,a,q,z,p=Ny-1,a-1,q=4,z,pTa,q,z→z′,pe-(Fy,a,q,z+Ma-1,q,p) when q=1Ny,a,q-1,z,pTa,q,z→z′,pe-(Fy,a,q,z+Ma,q,p) when q=2:4 6 SSBy,q,z,p=∑a=129Ny,a,q,z,pBa,pWa,p1000 7 Ny,a,q,s=∑West z=1:3East z=4:7Ny,a,q,z,p 8 Ny,a,q,z=∑West p=2East p=1Ny,a,q,z,p 9 Cy,a,s=(∑West z=1:3East z=4:7∑q=14∑p=12Ny,a,q,z,p Fy,a,q,zFy,a,q,z+Ma,q,p1-e-Fy,a,q,z+Ma,q,p)eεy,a,s 10 Iy,g,s=∑West z=1:3East z=4:7∑a=129∑p=12Sa,gNy,a,q,z,pWa,pQgeεy,g,s 11 Ey,g=∑aAxy,a,z,giy,g,s 12 Xy,a,g,s=(∑West z=1:3East z=4:7∑q=14∑p=12Ny,a,q,z,pEy,gQgSa,gEy,gQgSa,g+Ma,q,p1-e-Ey,gQgSa,g+Ma,q,p)eεy,a,s Estimation model 13 Ny,a=Fy,a+My,aFy,a1-e-Fy,a+My,aCy,a 14 Fy-1,a-1=Fy-1,a-1+My-1,a-1Cy-1,a-1Ny,aeFy-1,a-1+My-1,a-1-1 15 Ny,A=Fy,A+My,AFy,A1-e-Fy,A+My,ACy,A 16 Ny-1,A=Fy-1,A+My-1,AFy-1,A1-e-Fy-1,A+My-1,ACy-1,A 17 Fy-1,A-1=Fy-1,A-1+My-1,A-1Cy-1,A-1eFy-1,A-1+My-1,A-1-1Ny,A-Ny-1,AeFy-1,A+My-1,A 18 I^y,g,s=Qy,g,s∑aAΔy,a,g,sSy,a,g,sWy,a,g,sNy,a,s Performance metrics 19 relative bias= medianymediannθestn,y-θtrueyθtruey×100% 20 absolute bias=medianymediannθestn,y-θtruey Equation number . Operating model . 1 La,p=2=L1b+L2b-L1b1-e-Ka-A11-e-KA2-A11b 2 La,p=1=L∞1-e-Ka-a0 3 Wa=αLaβ 4 Fy,a,q,z=fy,a,sDz,qxy,a,z,g∑West z=1:3East z=4:7xy,a,z,g 5 Ny,a,q,z,p=Ny-1,a-1,q=4,z,pTa,q,z→z′,pe-(Fy,a,q,z+Ma-1,q,p) when q=1Ny,a,q-1,z,pTa,q,z→z′,pe-(Fy,a,q,z+Ma,q,p) when q=2:4 6 SSBy,q,z,p=∑a=129Ny,a,q,z,pBa,pWa,p1000 7 Ny,a,q,s=∑West z=1:3East z=4:7Ny,a,q,z,p 8 Ny,a,q,z=∑West p=2East p=1Ny,a,q,z,p 9 Cy,a,s=(∑West z=1:3East z=4:7∑q=14∑p=12Ny,a,q,z,p Fy,a,q,zFy,a,q,z+Ma,q,p1-e-Fy,a,q,z+Ma,q,p)eεy,a,s 10 Iy,g,s=∑West z=1:3East z=4:7∑a=129∑p=12Sa,gNy,a,q,z,pWa,pQgeεy,g,s 11 Ey,g=∑aAxy,a,z,giy,g,s 12 Xy,a,g,s=(∑West z=1:3East z=4:7∑q=14∑p=12Ny,a,q,z,pEy,gQgSa,gEy,gQgSa,g+Ma,q,p1-e-Ey,gQgSa,g+Ma,q,p)eεy,a,s Estimation model 13 Ny,a=Fy,a+My,aFy,a1-e-Fy,a+My,aCy,a 14 Fy-1,a-1=Fy-1,a-1+My-1,a-1Cy-1,a-1Ny,aeFy-1,a-1+My-1,a-1-1 15 Ny,A=Fy,A+My,AFy,A1-e-Fy,A+My,ACy,A 16 Ny-1,A=Fy-1,A+My-1,AFy-1,A1-e-Fy-1,A+My-1,ACy-1,A 17 Fy-1,A-1=Fy-1,A-1+My-1,A-1Cy-1,A-1eFy-1,A-1+My-1,A-1-1Ny,A-Ny-1,AeFy-1,A+My-1,A 18 I^y,g,s=Qy,g,s∑aAΔy,a,g,sSy,a,g,sWy,a,g,sNy,a,s Performance metrics 19 relative bias= medianymediannθestn,y-θtrueyθtruey×100% 20 absolute bias=medianymediannθestn,y-θtruey Equation elements are defined in Tables 1 and 3. Estimation model equations are adapted from Porch (2003). Open in new tab Table 3. Definitions and sources of equation elements used in Table 2. Element . Definition . Source . a Age n/a (subscript) p Population n/a (subscript) y Year n/a (subscript) q Quarter n/a (subscript) z Geographic zone n/a (subscript) s Stock n/a (subscript) g Fleet n/a (subscript) z→z' Movement from starting zone z to destination zone z' n/a (subscript) A Age (plus group) n/a (subscript) n Realization n/a (subscript) ^ Predicted (e.g. predicted index or catch) n/a (accent) θ Parameter value (e.g. recruitment, SSB, fishing mortality) n/a L Length (cm) Derived (Table 2 eq. 1 for the west and eq. 2 for the east) W Weight (kg) Derived (Table 2 eq. 3) F Fishing mortality rate (q−1 in operating model, y−1 in estimation model) Derived (Table 2 eq. 4) f Fishing mortality rate (y−1) Assumed (ICCAT, 2018a) D Quarterly proportion of fishing mortality rates by geographic zone Assumed (Supplementary Table S4) x Fleet-specific relative age-composition data by geographic zone Assumed (ICCAT, 2018a) N Abundance (numbers) Assumed (⁠ Ny=1974:2015,a=1,q=1,z=1,p=2 ⁠, Ny=1974:2015,a=1,q=1,z=7,p=1 ⁠, Ny=1974,a=1:29,q=1,z=1,p=2 ⁠, and Ny=1974,a=1:29,q=1,z=7,p=1 from Supplementary Tables S1 and S2); otherwise derived T Proportional movement of fish Assumed (Figure 2) M Natural mortality rate (q−1 in operating model, y−1 in estimation model) Assumed (Table 1) SSB Spawning stock biomass (tonnes) Derived (Table 2 eq. 6) B Spawning fraction Assumed (Table 1) C Catch (numbers) Derived (Table 2 eq. 9) ε Error Derived (Supplementary Table S5) I Index value Derived (Table 2 eq. 10 and 18) S Selectivity Assumed (Supplementary Tables S8 and S9) Q Catchability Assumed (Supplementary Table S10) E Effort Derived (Table 2 eq. 11 and Supplementary Tables S6 and S7) i Index value Assumed (ICCAT, 2018a) X Fleet-specific relative age-composition data by stock Derived (Table 2 eq. 12) σ Standard deviation Derived (Supplementary Table S5) Δ Adjustment for time of year (for predicted index values) Assumed (Tables 7 and 9 in ICCAT, 2018a) Element . Definition . Source . a Age n/a (subscript) p Population n/a (subscript) y Year n/a (subscript) q Quarter n/a (subscript) z Geographic zone n/a (subscript) s Stock n/a (subscript) g Fleet n/a (subscript) z→z' Movement from starting zone z to destination zone z' n/a (subscript) A Age (plus group) n/a (subscript) n Realization n/a (subscript) ^ Predicted (e.g. predicted index or catch) n/a (accent) θ Parameter value (e.g. recruitment, SSB, fishing mortality) n/a L Length (cm) Derived (Table 2 eq. 1 for the west and eq. 2 for the east) W Weight (kg) Derived (Table 2 eq. 3) F Fishing mortality rate (q−1 in operating model, y−1 in estimation model) Derived (Table 2 eq. 4) f Fishing mortality rate (y−1) Assumed (ICCAT, 2018a) D Quarterly proportion of fishing mortality rates by geographic zone Assumed (Supplementary Table S4) x Fleet-specific relative age-composition data by geographic zone Assumed (ICCAT, 2018a) N Abundance (numbers) Assumed (⁠ Ny=1974:2015,a=1,q=1,z=1,p=2 ⁠, Ny=1974:2015,a=1,q=1,z=7,p=1 ⁠, Ny=1974,a=1:29,q=1,z=1,p=2 ⁠, and Ny=1974,a=1:29,q=1,z=7,p=1 from Supplementary Tables S1 and S2); otherwise derived T Proportional movement of fish Assumed (Figure 2) M Natural mortality rate (q−1 in operating model, y−1 in estimation model) Assumed (Table 1) SSB Spawning stock biomass (tonnes) Derived (Table 2 eq. 6) B Spawning fraction Assumed (Table 1) C Catch (numbers) Derived (Table 2 eq. 9) ε Error Derived (Supplementary Table S5) I Index value Derived (Table 2 eq. 10 and 18) S Selectivity Assumed (Supplementary Tables S8 and S9) Q Catchability Assumed (Supplementary Table S10) E Effort Derived (Table 2 eq. 11 and Supplementary Tables S6 and S7) i Index value Assumed (ICCAT, 2018a) X Fleet-specific relative age-composition data by stock Derived (Table 2 eq. 12) σ Standard deviation Derived (Supplementary Table S5) Δ Adjustment for time of year (for predicted index values) Assumed (Tables 7 and 9 in ICCAT, 2018a) Open in new tab Table 3. Definitions and sources of equation elements used in Table 2. Element . Definition . Source . a Age n/a (subscript) p Population n/a (subscript) y Year n/a (subscript) q Quarter n/a (subscript) z Geographic zone n/a (subscript) s Stock n/a (subscript) g Fleet n/a (subscript) z→z' Movement from starting zone z to destination zone z' n/a (subscript) A Age (plus group) n/a (subscript) n Realization n/a (subscript) ^ Predicted (e.g. predicted index or catch) n/a (accent) θ Parameter value (e.g. recruitment, SSB, fishing mortality) n/a L Length (cm) Derived (Table 2 eq. 1 for the west and eq. 2 for the east) W Weight (kg) Derived (Table 2 eq. 3) F Fishing mortality rate (q−1 in operating model, y−1 in estimation model) Derived (Table 2 eq. 4) f Fishing mortality rate (y−1) Assumed (ICCAT, 2018a) D Quarterly proportion of fishing mortality rates by geographic zone Assumed (Supplementary Table S4) x Fleet-specific relative age-composition data by geographic zone Assumed (ICCAT, 2018a) N Abundance (numbers) Assumed (⁠ Ny=1974:2015,a=1,q=1,z=1,p=2 ⁠, Ny=1974:2015,a=1,q=1,z=7,p=1 ⁠, Ny=1974,a=1:29,q=1,z=1,p=2 ⁠, and Ny=1974,a=1:29,q=1,z=7,p=1 from Supplementary Tables S1 and S2); otherwise derived T Proportional movement of fish Assumed (Figure 2) M Natural mortality rate (q−1 in operating model, y−1 in estimation model) Assumed (Table 1) SSB Spawning stock biomass (tonnes) Derived (Table 2 eq. 6) B Spawning fraction Assumed (Table 1) C Catch (numbers) Derived (Table 2 eq. 9) ε Error Derived (Supplementary Table S5) I Index value Derived (Table 2 eq. 10 and 18) S Selectivity Assumed (Supplementary Tables S8 and S9) Q Catchability Assumed (Supplementary Table S10) E Effort Derived (Table 2 eq. 11 and Supplementary Tables S6 and S7) i Index value Assumed (ICCAT, 2018a) X Fleet-specific relative age-composition data by stock Derived (Table 2 eq. 12) σ Standard deviation Derived (Supplementary Table S5) Δ Adjustment for time of year (for predicted index values) Assumed (Tables 7 and 9 in ICCAT, 2018a) Element . Definition . Source . a Age n/a (subscript) p Population n/a (subscript) y Year n/a (subscript) q Quarter n/a (subscript) z Geographic zone n/a (subscript) s Stock n/a (subscript) g Fleet n/a (subscript) z→z' Movement from starting zone z to destination zone z' n/a (subscript) A Age (plus group) n/a (subscript) n Realization n/a (subscript) ^ Predicted (e.g. predicted index or catch) n/a (accent) θ Parameter value (e.g. recruitment, SSB, fishing mortality) n/a L Length (cm) Derived (Table 2 eq. 1 for the west and eq. 2 for the east) W Weight (kg) Derived (Table 2 eq. 3) F Fishing mortality rate (q−1 in operating model, y−1 in estimation model) Derived (Table 2 eq. 4) f Fishing mortality rate (y−1) Assumed (ICCAT, 2018a) D Quarterly proportion of fishing mortality rates by geographic zone Assumed (Supplementary Table S4) x Fleet-specific relative age-composition data by geographic zone Assumed (ICCAT, 2018a) N Abundance (numbers) Assumed (⁠ Ny=1974:2015,a=1,q=1,z=1,p=2 ⁠, Ny=1974:2015,a=1,q=1,z=7,p=1 ⁠, Ny=1974,a=1:29,q=1,z=1,p=2 ⁠, and Ny=1974,a=1:29,q=1,z=7,p=1 from Supplementary Tables S1 and S2); otherwise derived T Proportional movement of fish Assumed (Figure 2) M Natural mortality rate (q−1 in operating model, y−1 in estimation model) Assumed (Table 1) SSB Spawning stock biomass (tonnes) Derived (Table 2 eq. 6) B Spawning fraction Assumed (Table 1) C Catch (numbers) Derived (Table 2 eq. 9) ε Error Derived (Supplementary Table S5) I Index value Derived (Table 2 eq. 10 and 18) S Selectivity Assumed (Supplementary Tables S8 and S9) Q Catchability Assumed (Supplementary Table S10) E Effort Derived (Table 2 eq. 11 and Supplementary Tables S6 and S7) i Index value Assumed (ICCAT, 2018a) X Fleet-specific relative age-composition data by stock Derived (Table 2 eq. 12) σ Standard deviation Derived (Supplementary Table S5) Δ Adjustment for time of year (for predicted index values) Assumed (Tables 7 and 9 in ICCAT, 2018a) Open in new tab Operating model values for the resource size were represented both from a stock perspective, referring to the geographically distinct western and eastern mixed-population stocks separated by the 45°W stock boundary, and a population perspective, referring to the genetically distinct western and eastern populations originating from their respective natal grounds. Stock attributes were derived by summing over all geographic zones contained in a stock area (zones 1–3 for the western stock area and 4–7 for the eastern stock area; eq. 7; Tables 2 and 3). Population attributes were derived by summing over all fish that originated in the same spawning areas (Gulf of Mexico for the western population and Mediterranean Sea for the eastern population; eq. 8; Tables 2 and 3). The observation model generated pseudodata for catch-at-age, indices of relative abundance, and relative age composition of indices. Pseudodata for annual catch-at-age for each stock were derived using the Baranov catch equation (eq. 9; Tables 2 and 3). Pseudodata for fishery-dependent and fishery-independent indices of relative abundance were derived from simulated fishing fleets and surveys that emulated the geographic scope, magnitude, and time frame of index data available for the 2017 ICCAT stock assessments. These include catch-per-unit-effort of several commercial fleets for a range of gear types (rod and reel, longline, trap, baitboat), larval surveys in known spawning areas, aerial surveys, and an acoustic survey (ICCAT, 2017). Abundance values used to derive indices (eq. 10; Tables 2 and 3) were assumed to be from the beginning of the third quarter (to reflect the fall season when the majority of fishing effort occurs in the Atlantic bluefin tuna fishery). The only exception was for indices that measured relative abundance in spawning zones (zones 1 and 7), in which case the abundance was assumed to be from the beginning of the first quarter (to reflect the spawning season). The index-specific fishing mortality rate was derived as the product of the index-specific effort (eq. 11; Tables 2 and 3), catchability, and selectivity estimates from the 2017 ICCAT VPAs (Supplementary Tables S6–S10). This was used to derive the relative fishing mortality rate-at-age within fleets and years in the calculation of pseudodata for the relative age composition of indices (eq. 12; Tables 2 and 3). The observation model assumed age plus groups of 1–16+ for the western stock and 1–10+ for the eastern stock, as in ICCAT (2018a; Table 1). Observation error in log space, where ε∼N(0;σ2) ⁠, was generated for the catch-at-age of each stock, the relative age composition of each index, and each index of relative abundance. Observation error for the catch-at-age and relative age composition, εy,a,s ⁠, had a stock-specific standard deviation σs calculated as the root mean square error of the observed and predicted catches from an exploratory age-structured assessment programme (Legault and Restrepo, 1998) analysis of Atlantic bluefin tuna data (Maguire et al., 2018; Supplementary Table S5). Observation error for each index, εy,g,s ⁠, had a standard deviation σg calculated as the root mean square error of the observed and predicted index values from the 2017 ICCAT western and eastern bluefin tuna VPAs (ICCAT, 2018a; Supplementary Table S5). Estimation models We tested the performance of the VPA-2BOX programme (version 4.01; Porch et al., 2001) for estimating resource size, fishing mortality rates, and F0.1 reference points. VPA (as opposed to other age-structured models such as statistical catch-at-age) was used in this study because it has been the status quo for Atlantic bluefin tuna stock assessment for over two decades. VPA-2BOX is a calibrated VPA, which is an age-structured model that estimates terminal fishing mortality rates (terminal F) by fitting indices to recursive abundance calculations derived from terminal F. However, similar to more integrated analyses, VPA-2BOX also includes fleet structure to inform selectivity of fishery catch rate indices (Porch, 2003). Model settings were based on ICCAT (2018a) western and eastern Atlantic bluefin tuna stock assessment base case runs. Selected equations are documented in Table 2 for the backwards recursion calculating abundance and fishing mortality for all non-plus group ages (eq. 13, 14; Tables 2 and 3) and plus group ages (eq. 15-17), and for calculating predicted indices (eq. 18; see Porch, 2003 for further details). The results of the western and eastern estimation models were used to calculate F0.1 reference points, a proxy for FMSY ⁠. To specifically test the current ICCAT management process (ICCAT Standing Committee on Research and Statistics (SCRS), 2017), estimated values of F0.1 were compared to the estimated current fishing mortality rate, Fcurrent ⁠, to determine current exploitation status. If Fcurrent/F0.1 was >1, overfishing was occurring; if ≤1, overfishing was not occurring (see Supplementary data B for details on calculation of F0.1) ⁠. Performance evaluation We evaluated the performance of the closed-population western and eastern VPAs on pseudodata generated with fish movement and population mixing. To evaluate the mismatch in operating model and estimation model structural assumptions, we conducted “cross-tests” and “self-tests” (Deroba et al., 2015). Cross-tests reveal changes in estimation performance due to structural differences between the operating model and estimation model. Cross-tests were performed by applying the VPA estimation model, which assumes no spatial structure and annual time steps, to pseudodata generated by the spatially complex, seasonally structured operating model. To isolate the effects of spatial misspecification, population parameter settings other than spatial and temporal structure were the same between the operating model and estimation model, including maturity, natural mortality, and weight-at-age. In contrast, self-tests require an estimation model that makes similar structural assumptions to the operating model and facilitate understanding of sources of additional estimation bias other than structural misspecification. In our self-tests, the operating model was simplified (fish movement eliminated, seven geographic zones aggregated into two stock areas, quarterly time steps converted to annual time steps) to generate pseudodata reflective of the single-stock, single-area structural assumptions of the VPA model (i.e. producing stock-level biomass only, not population-level). For each cross-test and self-test of western and eastern fish, 500 realizations were run and estimation model results compared to the true values from the operating model. Relative and absolute bias of median estimated spawning stock biomass (SSB), recruitment, and F/F0.1 across all 42 years of all 500 simulated realizations were calculated relative to the true values from the operating model (eq. 19, 20; Tables 2 and 3). Bias in estimated Fcurrent/F0.1 was calculated once at the end of each time series across all realizations test the current ICCAT management process. In both cross-tests and self-tests, bias was calculated relative to the true stock values from the operating model. Because operating model spatial structure in the cross-tests produced different population and stock perspectives, in cross-tests only, bias was also calculated for estimated SSB, F/F0.1 ⁠, and Fcurrent/F0.1 relative to the true population values from the operating model to reflect the closed-population assumptions of the VPA. Results Operating model dynamics Fish movement resulted in a divergence between population and stock SSB over time after the 1974 model initialization (Figure 3c and d). There was net movement of eastern population fish into the western stock area, resulting in a larger western stock than western population, and a smaller eastern stock than eastern population (Figure 3e and f). The difference in size between the eastern stock and population was slight, but the western stock was larger than the western population by as much as 300% in some years. SSB and recruitment of the eastern population were greater than the western population throughout the time series by approximately an order of magnitude (Figure 3a–d). Despite the eastern population being the largest contribution to SSB in both the western and eastern stock areas, the majority of western and eastern population SSB was found within their respective stock areas for both populations (Figure 3e and f). Figure 3. Open in new tabDownload slide Operating model time series of western and eastern population recruitment (a and b), population and stock SSB (c and d), and population SSB in the western and eastern stock areas (e and f) for the base case movement scenario. SSB values were taken from the fall season (quarter = 3). Figure 3. Open in new tabDownload slide Operating model time series of western and eastern population recruitment (a and b), population and stock SSB (c and d), and population SSB in the western and eastern stock areas (e and f) for the base case movement scenario. SSB values were taken from the fall season (quarter = 3). Fishing mortality rates were generally greater, and abundance was generally lower in the western stock area relative to the eastern stock area (Supplementary Figures S4 and S6). The realized fishing mortality rates in the operating model were lower than the ICCAT VPA results as a result of conditioning on the 2017 ICCAT VPAs, partitioning fishing mortality to season and geographic zones according to fishing patterns, and seasonally distributing fish to areas (Supplementary Figure S4). The lower effective fishing mortality rates resulted in higher SSB in the operating model compared to the 2017 ICCAT VPA estimates. Current efforts by ICCAT scientists to condition operating models for management strategy evaluation (MSE) of Atlantic bluefin tuna have had similar difficulties reflecting the ICCAT stock assessment perception of stock development when incorporating stock mixing (Carruthers and Kell, 2017; Carruthers and Butterworth, 2018a, b). See Supplementary data C for further details on operating model dynamics. True values of F/F0.1 for the entire time series and Fcurrent/F0.1 at the end of the time series were nearly identical between the operating models for the self-test and the stock perspective of the base case movement cross-test (Figures 4e and f and 5e and f). True F/F0.1 values for the self-test and cross-test were within 0.04 of each other for the eastern stock, and within 0.005 of each other for the western stock across the entire time series. True Fcurrent/F0.1 values for the self-test and cross-test were within 0.002 of each other for both stocks. In the cross-test operating model, Fcurrent/F0.1 was relatively higher for the western stock (0.45) than for the eastern stock (0.12). In contrast, the true population values for these quantities were lower (0.23 for the western population and 0.05 for the eastern population), indicating a very low probability of overfishing. Figure 4. Open in new tabDownload slide Results of the self-tests for recruitment (a and b), SSB (c and d), and F/F0.1 (e and f) by stock. Boxplots display median, interquartile range (IQR, hinges), and 1.5 × IQR (whiskers) for the estimation model results. Solid black lines are the operating model, and solid horizontal black lines in (e) and (f) are the F/F0.1=1 line (above which overfishing is occurring). Outliers and poorly estimated recruitment values from 2012 to 2015 are not shown. Figure 4. Open in new tabDownload slide Results of the self-tests for recruitment (a and b), SSB (c and d), and F/F0.1 (e and f) by stock. Boxplots display median, interquartile range (IQR, hinges), and 1.5 × IQR (whiskers) for the estimation model results. Solid black lines are the operating model, and solid horizontal black lines in (e) and (f) are the F/F0.1=1 line (above which overfishing is occurring). Outliers and poorly estimated recruitment values from 2012 to 2015 are not shown. Performance of estimation models West Median bias in western recruitment increased from 22% in the self-test (Table 4 and Figure 4a) to 147% in the cross-test (Table 4 and Figure 5a), but the estimation model replicated some relative trends in recruitment from the operating model. SSB was more accurately estimated in the western self-test with only −10% median bias (Table 4 and Figure 4c), but median bias increased in the cross-test to 92% relative to the true population SSB and to −48% relative to the true stock SSB (Table 4 and Figure 5c). SSB estimates from the cross-test mostly fell between the true operating model values for the western population and stock SSB (Figure 5c). The median bias in F/F0.1 was 23% in the western self-test (Table 4 and Figure 4e), but increased slightly in the cross-test to 57% relative to the true population F/F0.1 and decreased to −13% relative to the true stock F/F0.1 (Table 4 and Figure 5e). The median bias in Fcurrent/F0.1 was 49% in the western self-test, but increased in the cross-test to 113% relative to the true population Fcurrent/F0.1 and decreased to 10% relative to the true stock Fcurrent/F0.1 (Table 4). Figure 5. Open in new tabDownload slide Results of the cross-tests for recruitment (a and b), SSB (c and d), and F/F0.1 (e and f) by population/stock. Boxplots display median, IQR (hinges), and 1.5 × IQR (whiskers) for the estimation model results. Solid black lines are the operating model population, dashed black lines are the operating model stock, and solid horizontal black lines in (e) and (f) are the F/F0.1=1 line (above which overfishing is occurring). Outliers and poorly estimated recruitment values from 2012 to 2015 are not shown. Figure 5. Open in new tabDownload slide Results of the cross-tests for recruitment (a and b), SSB (c and d), and F/F0.1 (e and f) by population/stock. Boxplots display median, IQR (hinges), and 1.5 × IQR (whiskers) for the estimation model results. Solid black lines are the operating model population, dashed black lines are the operating model stock, and solid horizontal black lines in (e) and (f) are the F/F0.1=1 line (above which overfishing is occurring). Outliers and poorly estimated recruitment values from 2012 to 2015 are not shown. Table 4. Relative and absolute median bias (in parentheses) in VPA model estimates of recruitment (1974–2011), SSB (1974–2015), F/F0.1 (1976–2015), and Fcurrent/F0.1 (2012–2014) relative to true values from the operating model for the self-tests and cross-tests. Unit . Test . Group . Recruitment (numbers) . SSB (tonnes) . F/F0.1 . Fcurrent/F0.1 . Self-test Stock 22% (34 270) −10% (−1 132) 23% (0.26) 49% (0.22) West Cross-test Population 147% (230 960) 92% (19 597) 57% (0.31) 113% (0.26) Stock −48% (−32 821) −13% (−0.15) 10% (0.04) Self-test Stock 15% (254 626) −44% (−171 186) 54% (0.31) 171% (0.21) East Cross-test Population −46% (−1 126 789) −81% (−507 751) 204% (0.59) 481% (0.25) Stock −80% (−461 161) 23% (0.16) 152% (0.18) Unit . Test . Group . Recruitment (numbers) . SSB (tonnes) . F/F0.1 . Fcurrent/F0.1 . Self-test Stock 22% (34 270) −10% (−1 132) 23% (0.26) 49% (0.22) West Cross-test Population 147% (230 960) 92% (19 597) 57% (0.31) 113% (0.26) Stock −48% (−32 821) −13% (−0.15) 10% (0.04) Self-test Stock 15% (254 626) −44% (−171 186) 54% (0.31) 171% (0.21) East Cross-test Population −46% (−1 126 789) −81% (−507 751) 204% (0.59) 481% (0.25) Stock −80% (−461 161) 23% (0.16) 152% (0.18) Open in new tab Table 4. Relative and absolute median bias (in parentheses) in VPA model estimates of recruitment (1974–2011), SSB (1974–2015), F/F0.1 (1976–2015), and Fcurrent/F0.1 (2012–2014) relative to true values from the operating model for the self-tests and cross-tests. Unit . Test . Group . Recruitment (numbers) . SSB (tonnes) . F/F0.1 . Fcurrent/F0.1 . Self-test Stock 22% (34 270) −10% (−1 132) 23% (0.26) 49% (0.22) West Cross-test Population 147% (230 960) 92% (19 597) 57% (0.31) 113% (0.26) Stock −48% (−32 821) −13% (−0.15) 10% (0.04) Self-test Stock 15% (254 626) −44% (−171 186) 54% (0.31) 171% (0.21) East Cross-test Population −46% (−1 126 789) −81% (−507 751) 204% (0.59) 481% (0.25) Stock −80% (−461 161) 23% (0.16) 152% (0.18) Unit . Test . Group . Recruitment (numbers) . SSB (tonnes) . F/F0.1 . Fcurrent/F0.1 . Self-test Stock 22% (34 270) −10% (−1 132) 23% (0.26) 49% (0.22) West Cross-test Population 147% (230 960) 92% (19 597) 57% (0.31) 113% (0.26) Stock −48% (−32 821) −13% (−0.15) 10% (0.04) Self-test Stock 15% (254 626) −44% (−171 186) 54% (0.31) 171% (0.21) East Cross-test Population −46% (−1 126 789) −81% (−507 751) 204% (0.59) 481% (0.25) Stock −80% (−461 161) 23% (0.16) 152% (0.18) Open in new tab When mean and low movement rates were compared in the cross-test, estimates of western recruitment, population SSB, and stock SSB in the low movement scenario were less biased (32, −12, and −43%, respectively; Supplementary Figure S3a and c) than in the base case movement scenario (147, 92, and −48% respectively; Supplementary Table S11). Estimates of population and stock F/F0.1 were similarly biased in the low movement scenario (68 and −17%, respectively; Supplementary Figure S3e) as in the base case scenario (57 and −13%, respectively; Supplementary Table S11). Estimates of population and stock Fcurrent/F0.1 were similarly biased in the low movement scenario (130 and 10%, respectively) as in the base case scenario (113 and 10%, respectively; Supplementary Table S11). East Median bias in eastern recruitment was 15% in the self-test (Table 4 and Figure 4b) and −46% in the cross-test (Table 4 and Figure 5b). The median bias in eastern SSB was −44% in the self-test (Table 4 and Figure 4d), and increased in the cross-test to −81 and −80% relative to the true population and stock SSB, respectively (Table 4 and Figure 5d). Median bias in F/F0.1 was 54% in the eastern self-test (Table 4 and Figure 4f) but increased in the cross-test to 204% relative to the true population F/F0.1 and decreased to 23% relative to the true stock F/F0.1 (Table 4 and Figure 5f). The median bias in Fcurrent/F0.1 was 171% in the eastern self-test but increased in the cross-test to 481% relative to the true population Fcurrent/F0.1 and decreased to 152% relative to the true stock Fcurrent/F0.1 (Table 4). When mean and low movement rates were compared in the cross-test, estimates of eastern recruitment, population SSB, and stock SSB in the low movement scenario were slightly more biased (−60, −85, and −85%, respectively; Supplementary Figure S3b and d) than in the base case movement scenario (−46, −81, and −80%, respectively; Supplementary Table S11). Estimates of population and stock F/F0.1 were similarly biased in the low movement scenario (191 and 15%, respectively; Supplementary Figure S3f) as in the base case scenario (204 and 23%, respectively; Supplementary Table S11). Estimates of eastern population and stock Fcurrent/F0.1 were strongly biased in both the low movement scenario (685 and 194%, respectively) and the base case scenario (481 and 152%, respectively; Supplementary Table S11). Discussion Correctly specifying models is a challenge to all stock assessment scientists, especially for fisheries with complex spatial structure (Berger et al., 2017; Punt, 2019a). For species like Atlantic bluefin tuna, in which two distinct populations have strong connectivity, a single-stock, single-area stock assessment model fit to data from each mixed-population stock is inherently misspecified (Goethel et al., 2015; Jardim et al., 2018). Non-spatial assessments that aggregate spatially structured data have been shown by this study and others to result in biased estimates of stock size (e.g. Cope and Punt, 2011; Ying et al., 2011; Punt et al., 2015). Stock assessment approaches that can account for (or are robust to) the dynamics of movement and mixing are needed. In this study, recruitment and SSB of western Atlantic bluefin tuna were overestimated due to the substantial subsidies of the much larger eastern population into the western stock area. Management procedures based on absolute estimates of the western stock size derived from single-stock, single-area models could pose a high risk of jeopardizing the sustainability of the smaller western population (Kerr et al., 2017, 2018). On the other hand, eastern SSB was underestimated relative to both the true population and stock, implying that management procedures based on absolute stock size estimates for the larger eastern resource may be robust to stock mixing and effective at preventing overfishing. However, there is trade-off in the loss of potential fishery yield. Vincent et al. (2017) also found that stock assessment models may have difficulty estimating individual population abundances under certain movement rates when productivities of intermixing populations were vastly different. Earlier studies (for salmon, Ricker, 1958; and Atlantic cod, Fu and Fanning, 2004) have demonstrated the unintended consequences of managing mixed stocks without consideration for the distinct population dynamics and productivities that make up these stocks. Even if catch limits are sustainable from the stock perspective of intermixed populations, there is a risk of overfishing the smaller population because its signal may be washed out by the larger population, resulting in depletion or extermination of the smaller population. In addition to structural model misspecification, this VPA model application has a demonstrated history of poor performance and instability that has been attributed to poor data quality (ICCAT, 2015, 2018a) and unsuitable parameter settings (Zarrad et al., 2018) for Atlantic bluefin tuna, which are common limitations to fisheries stock assessment (Magnusson and Hilborn, 2007). These concerns are magnified in a simulation scenario where human aspects of the stock assessment modelling process, such as iterative tuning of parameter starting values to fit the model to the data, cannot be automated for each realization. Because of these existing concerns with the VPA model, we conducted the self-test in contrast to the base case movement cross-test to identify whether estimation bias could have resulted from the instability and poor performance of the VPA instead of (or in addition to) structural model misspecification. Relative to the cross-tests, both the western and eastern self-tests produced estimates of recruitment and SSB with low bias, confirming that structural model misspecification did account for most of the bias in these estimates in the cross-tests. In addition, comparisons between the base case movement and low movement cross-tests indicated whether the magnitude of movement influenced estimation performance. Estimates of eastern recruitment and SSB did not seem to become more accurate at lower movement rates, which suggest that the presence, not the magnitude, of movement caused bias in these estimates. On the other hand, western recruitment and population SSB were more accurately estimated when movement rates were lower, because the smaller western population was less washed out by the larger eastern population. For all operating model scenarios tested in this study, we found that the ICCAT management reference value of Fcurrent/F0.1 nearly always identified the correct exploitation status (overfishing not occurring), under the current conditions for Atlantic bluefin tuna. Current conditions include the perception of low exploitation as modelled after the 2017 ICCAT stock assessments. However, this metric overestimated the true value most of the time and particularly strongly relative to the true population values. This was generally also true for F/F0.1 for all preceding years of the time series (except for its being underestimated relative to the true western stock values). Because of the strong estimation bias in this metric, it might only be considered reliable as a reference point when Fcurrent/F0.1≪1 or ≫1 ⁠. Conditioning the operating model in a simulation study so that it represents the system of interest is particularly important. Previous efforts to fit spatially explicit tag-integrated stock assessment models to Atlantic bluefin tuna data [MAST model documented in Taylor et al. (2011) and Kerr et al. (2017)] were not stable enough to provide the flexibility or stability to support the alternative scenarios needed for this evaluation (Kerr et al., 2018). Instead, we conditioned the operating model on recent stock assessments (ICCAT, 2018a) and simulated mixing based on telemetry (Galuardi et al., 2014). The operating model framework documented here built upon preceding studies (Kerr et al., 2017, 2018) supports alternative life history scenarios, which facilitates conditioning of future simulation studies to reflect uncertainties in parameter assumptions (e.g. natural mortality, spawning fraction; e.g. Weston et al., 2019). Other spatial simulation studies for Atlantic cod (Kerr et al., 2014) and Pacific bluefin tuna (Lee et al., 2017) have effectively conditioned operating models on results from simpler estimation models with modifications to reflect the plausible spatial structure. Alternative approaches have attempted to overcome the limitations of using stock-level data for informing sustainable management of mixed populations. For example, attempts to estimate population-of-origin by using stock composition data derived from otolith chemistry analyses to revise model input data produced similar stock assessment results for Atlantic bluefin tuna (Morse et al., 2018; see a similar approach in Li et al., 2018 for lake whitefish). In this study we modelled “white noise” observation error, in which the error terms were generated independently with constant variance, but alternative approaches to error generation could be considered. For all pseudodata types, variance was based on the observed error in Atlantic bluefin tuna catch and index data and was constant across ages and years. This approach allowed us to focus on the effect of structural error and to remain consistent with ICCAT approaches to observation error in MSE (Butterworth et al., 2016). However, the assumption of constant variance across ages and years, though computationally convenient, is not necessarily realistic. For example, when certain age classes are targeted more heavily than others, sample sizes of under-targeted age classes can be smaller and therefore observation errors higher. In future applications of this research (e.g. MSE) it would be valuable to test alternative scenarios of more “realistic” error structure in which the error is autocorrelated with age and time, i.e. “red noise” (e.g. Wiedenmann et al., 2015). This study demonstrated the limitations of using single-stock assessment models on mixed-population stocks, where the smaller population is at risk of overexploitation and the large population at risk of foregone yield. We hope this approach for addressing the complexities of spatial dynamics using a spatial simulation model may be applied to other fisheries where mixing of distinct population units has been identified but not operationalized in management (see case studies reviewed by Berger et al., 2017; Cadrin et al., 2019; Punt, 2019b). Another important consideration in addition to incorrectly specified model structure is incomplete, unreliable, or conflicting data sets, which have been cited as challenges to the practicality of spatially structured modelling approaches for mixed-population stocks (e.g. Thorrold et al., 2001; Goethel et al., 2015). Mixed-population stocks can be managed at a finer scale, in real time, and with spatio-temporal complexity (Dunn et al., 2011) with the development of more reliable data on mixing dynamics. Methods that can provide information on mixing dynamics include large-scale fishery-independent and fishery-dependent tagging studies, and stock composition data from otolith chemistry and genetic markers. These advancements could particularly aid efforts to monitor and protect the smaller, more vulnerable population in a mixed stock. Our findings justify further development of stock assessment approaches that explicitly account for stock mixing to reduce estimation bias resulting from model misspecification. Our findings elucidate limitations of the status quo stock assessment approach for Atlantic bluefin tuna, in which misspecified spatial structure creates inaccurate perceptions of population sizes and productivities. Different models that account for stock mixing have been considered by the ICCAT Bluefin Tuna Working Group but are still in development (Carruthers et al., 2016; ICCAT, 2018a; Morse et al., 2018). Next steps for research and management include the use of MSE (Punt et al., 2016) to test the long-term performance of different management procedures on mixed-population stocks. Extension of this research using MSE could test both model-based and empirical-based management procedures, as well as additional dimensions of uncertainty in population dynamics and life history, such as natural mortality, stock-recruitment relationships, and time-varying movement (e.g. Carruthers and Butterworth, 2018a, b; ICCAT, 2018b; Kerr et al., 2018; Weston et al., 2019). Acknowledgements We acknowledge Gavin Fay for sharing his expertise and insight to improve the simulation methodology and rigour of this study. We thank the following for contributing tagging data that informed movement parameterization: Large Pelagics Research Center, Fisheries and Oceans Canada, ICCAT Atlantic-Wide Research Program for Bluefin Tuna (GBYP), AZTI Tecnalia, and NOAA. We are grateful to Matt Lauretta, John Walter, Lee Qi, and the anonymous reviewers who provided feedback to improve this article. The many scientists of the ICCAT Bluefin Tuna Working Group contributed to this research, especially Clay Porch, Tristan Rouyer, Ai Kimoto, J.-J. Maguire, Dave Secor, and Doug Butterworth. Funding This research was funded by the United States National Oceanic and Atmospheric Administration (NOAA) Bluefin Tuna Research Program grant award NA16NMF4720101. The scientific results and conclusions, as well as any views or opinions expressed herein, are those of the authors and do not necessarily reflect those of NOAA or the Department of Commerce. References Berger A. M. , Goethel D. R., Lynch P. D., Quinn T. D., Mormede S., McKenzie J., Dunn A. 2017 . Space oddity: the mission for spatial integration . Canadian Journal of Fisheries and Aquatic Sciences , 74 : 1698 – 1716 . Google Scholar Crossref Search ADS WorldCat Block B. A. , Teo S. L. H., Walli A., Boustany A., Stokesbury M. J. W., Farwell C. J., Weng K. C., et al. 2005 . Electronic tagging and population structure of Atlantic bluefin tuna . Nature , 434 : 1121 – 1127 . Google Scholar Crossref Search ADS PubMed WorldCat Boustany A. M. , Reeb C. A., Block B. A. 2008 . Mitochondrial DNA and electronic tagging reveal population structure of Atlantic bluefin tuna (Thunnus thynnus) . Marine Biology , 156 : 13 – 24 . Google Scholar Crossref Search ADS WorldCat Butterworth D. S. , Carruthers T. R., Rademeyer R. A. 2016 . Specifications for MSE trials for Bluefin tuna in the North Atlantic. Document prepared for the ICCAT Atlantic Bluefin MSE meeting in Monterey, USA. Butterworth D. S. , Punt A. E. 1994 . The robustness of estimates of stock status for the western North Atlantic bluefin tuna population to violations of the assumptions underlying the associated assessment models . ICCAT Collective Volume of Scientific Papers , 42 : 192 – 210 . OpenURL Placeholder Text WorldCat Cadrin S. X. , Secor D. H. 2009 . Accounting for spatial population structure in stock assessment: Past, present, and future. In The Future of Fisheries Science in North America , pp. 405 – 426 . Ed. by Beamish R. J., Rothschild B. J.. Springer Science + Business Media B.V ., Dordrecht. 736 pp. Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Cadrin S. X. , Goethel D. R., Morse M. R., Fay G., Kerr L. A. 2019 . “So, where do you come from?” The impact of assumed spatial population structure on estimates of recruitment . Fisheries Research , 217 : 156 – 168 . Google Scholar Crossref Search ADS WorldCat Carruthers T. , Butterworth D. 2018 a. ABT-MSE: an R package for Atlantic bluefin tuna management strategy evaluation . ICCAT Collective Volume of Scientific Papers , 74 : 3553 – 3559 . OpenURL Placeholder Text WorldCat Carruthers T. , Butterworth D. 2018 b. Summary of a reference set of conditioned operating models for Atlantic bluefin tuna . ICCAT Collective Volume of Scientific Papers , 73 : 3523 – 3541 . OpenURL Placeholder Text WorldCat Carruthers T. , Kell L. 2017 . Issues arising from the preliminary conditioning of operating models for Atlantic bluefin tuna . ICCAT Collective Volume of Scientific Papers , 73 : 2171 – 2182 . OpenURL Placeholder Text WorldCat Carruthers T. , Kimoto A., Powers J. E., Kell L., Butterworth D. S., Lauretta M. V., Kitakado T. 2016 . Structure and estimation framework for Atlantic bluefin tuna operating models . ICCAT Collective Volume of Scientific Papers , 72 : 1782 – 1795 . OpenURL Placeholder Text WorldCat Cope J. M. , Punt A. E. 2011 . Reconciling stock assessment and management scales under conditions of spatially varying catch histories . Fisheries Research , 107 : 22 – 38 . Google Scholar Crossref Search ADS WorldCat Deroba J. J. , Butterworth D. S., Methot R. D. Jr, De Oliveira J. A. A., Fernandez C., Nielsen A., Cadrin S. X., et al. 2015 . Simulation testing the robustness of stock assessment models to error: some results from the ICES strategic initiative on stock assessment methods . ICES Journal of Marine Science , 72 : 19 – 30 . Google Scholar Crossref Search ADS WorldCat Dunn D. C. , Boustany A. M., Halpin P. N. 2011 . Spatio-temporal management of fisheries to reduce by-catch and increase fishing selectivity . Fish and Fisheries , 12 : 110 – 119 . Google Scholar Crossref Search ADS WorldCat Fromentin J.-M. 2006 . ICCAT Manual Chapter 2.1.5: Atlantic Bluefin Tuna. ICCAT, Madrid, Spain. https://www.iccat.int/Documents/SCRS/Manual/CH2/2_1_5_BFT_ENG.pdf (last accessed 1 March 2018). Fromentin J.-M. , Powers J. E. 2005 . Atlantic bluefin tuna: population dynamics, ecology, fisheries and management . Fish and Fisheries , 6 : 281 – 306 . Google Scholar Crossref Search ADS WorldCat Fu C. , Fanning P. 2004 . Spatial considerations in the management of Atlantic cod off Nova Scotia, Canada . North American Journal of Fisheries Management , 24 : 775 – 784 . Google Scholar Crossref Search ADS WorldCat Galuardi B. , Cadrin S. X., Arregi I., Arrizabalaga H., DiNatale A., Brown C., Lauretta M., et al. 2018 . Atlantic bluefin tuna area transition matrices estimated from electronic tagging and SATTAGSIM . ICCAT Collective Volume of Scientific Papers , 74 : 2903 – 2921 . OpenURL Placeholder Text WorldCat Galuardi B. , Cadrin S. X., Kerr L., Miller T. J., Lutcavage M. 2014 . Using electronic tag data to provide transition matrices for movement inclusive population models . ICCAT Collective Volume of Scientific Papers , 71 : 1746 – 1756 . OpenURL Placeholder Text WorldCat Goethel D. R. , Kerr L. A., Cadrin S. X. 2016 . Incorporating spatial population structure into the assessment-management interface of marine resources. In Management Science in Fisheries: An Introduction to Simulation-Based Methods , pp. 319 – 347 . Ed. by C. T. T Edwards, Dankel D. J.. Routledge , New York . 460 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Goethel D. R. , Legault C. M., Cadrin S. X. 2015 . Testing the performance of a spatially explicit tag-integrated stock assessment model of yellowtail flounder (Limanda ferruginea) through simulation analysis . Canadian Journal of Fisheries and Aquatic Sciences , 72 : 582 – 601 . Google Scholar Crossref Search ADS WorldCat Guan W. , Cao J., Chen Y., Cieri M. 2013 . Impacts of population and fishery spatial structures on fishery stock assessment . Canadian Journal of Fisheries and Aquatic Sciences , 70 : 1178 – 1189 . Google Scholar Crossref Search ADS WorldCat ICCAT. 1981 . Recommendation on Bluefin Management Measures 81-01. ICCAT, Madrid, Spain. https://www.iccat.int/Documents/Recs/compendiopdf-e/1981-01-e.pdf (last accessed 1 November 2017). ICCAT. 2015 . Report of the 2014 Atlantic bluefin tuna stock assessment session . ICCAT Collective Volume of Scientific Papers , 71 : 692 – 945 . OpenURL Placeholder Text WorldCat ICCAT. 2017 . Report of the 2017 ICCAT Bluefin Data Preparatory Meeting. ICCAT, Madrid, Spain. ICCAT. 2018 a. Report of the 2017 ICCAT bluefin stock assessment meeting . ICCAT Collective Volume of Scientific Papers , 74 : 2372 – 2535 . OpenURL Placeholder Text WorldCat ICCAT. 2018 b. Report of the 2018 ICCAT Bluefin Tuna Species Group MSE Intersessional Meeting. ICCAT, Madrid, Spain. ICCAT Standing Committee on Research and Statistics (SCRS). 2017 . Report of the Standing Committee on Research and Statistics (SCRS). ICCAT, Madrid, Spain. Jardim E. , Eero M., Silva A., Ulrich C., Pawlowski L., Holmes S. J., Ibaibarriaga L., et al. 2018 . Testing spatial heterogeneity with stock assessment models . PLoS ONE , 13 : e0190791 . Google Scholar Crossref Search ADS PubMed WorldCat Kerr L. A. , Cadrin S. X., Kovach A. I. 2014 . Consequences of a mismatch between biological and management units on our perception of Atlantic cod off New England . ICES Journal of Marine Science , 71 : 1366 – 1381 . Google Scholar Crossref Search ADS WorldCat Kerr L. A. , Cadrin S. X., Secor D. H., Taylor N. 2017 . Modeling the implications of stock mixing and life history uncertainty of Atlantic bluefin tuna . Canadian Journal of Fisheries and Aquatic Sciences , 74 : 1990 – 2004 . Google Scholar Crossref Search ADS WorldCat Kerr L. A. , Goethel D. R. 2014 . Simulation modeling as a tool for synthesis of stock identification information. In Stock Identification Methods: Applications in Fishery Science , 2nd edn, pp. 502 – 533 . Ed. by Cadrin S. X., Kerr L. A., Mariani S.. Elsevier , Waltham . 588 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Kerr L. A. , Morse M. R., Cadrin S. X., Galuardi B. 2018 . Application of an Atlantic bluefin tuna operating model to generate pseudodata for stock assessment testing . ICCAT Collective Volume of Scientific Papers , 74 : 3405 – 3426 . OpenURL Placeholder Text WorldCat Legault C. M. , Restrepo V. R. 1998 . A flexible forward age-structured assessment program . ICCAT Collective Volume of Scientific Papers , 49 : 246 – 253 . OpenURL Placeholder Text WorldCat Lee H. , Piner K. R., Maunder M. R., Taylor I. G., Methot R. D. Jr., 2017 . Evaluation of alternative modelling approaches to account for spatial effects due to age-based movement . Canadian Journal of Fisheries and Aquatic Sciences , 74 : 1832 – 1844 . Google Scholar Crossref Search ADS WorldCat Li Y. , Bence J. R., Brenden T. O. 2018 . Can spawning origin information of catch or a recruitment penalty improve assessment and fishery management performance for a spatially structured stock assessment model? Canadian Journal of Fisheries and Aquatic Sciences , 75 : 2136 – 2148 . Google Scholar Crossref Search ADS WorldCat Magnusson A. , Hilborn R. 2007 . What makes fisheries data informative? Fish and Fisheries , 8 : 337 – 358 . Google Scholar Crossref Search ADS WorldCat Maguire J.-J. , Cadrin S. X., Hanke A., Melvin G. 2018 . An exploration of bluefin tuna data in the North Atlantic (West and East + Mediterranean) with ASAP . ICCAT Collective Volume of Scientific Papers , 74 : 3205 – 3223 . OpenURL Placeholder Text WorldCat Morse M. , Cadrin S., Kerr L., Secor D., Siskey M., Arrizabalaga H., Hanke A., et al. 2018 . An updated analysis of bluefin tuna stock mixing . ICCAT Collective Volume of Scientific Papers , 74 : 3486 – 3509 . OpenURL Placeholder Text WorldCat National Research Council (NRC). 1994 . An Assessment of Atlantic Bluefin Tuna . National Academy Press , Washington, D.C . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Nemerson D. , Berkeley S., Safina C. 2000 . Spawning site fidelity in Atlantic bluefin tuna, Thunnus thynnus: the use of size-frequency analysis to test for the presence of migrant east Atlantic bluefin tuna on Gulf of Mexico spawning grounds . Fishery Bulletin , 98 : 118 – 126 . OpenURL Placeholder Text WorldCat Pons M. , Melnychuk M. C., Hilborn R. 2018 . Management effectiveness of large pelagic fisheries in the high seas . Fish and Fisheries , 19 : 260 – 270 . Google Scholar Crossref Search ADS WorldCat Porch C. E. 1995 . A two-area VPA with discrete mixing: can we discriminate between mixing rates given the present condition of the data? ICCAT Collective Volume of Scientific Papers , 44 : 198 – 208 . OpenURL Placeholder Text WorldCat Porch C. E. 2003 . VPA-2BOX Version 3.01 User’s Guide. NOAA Sustainable Fisheries Division Contribution SFD/2003-2004, Miami, FL. Porch C. E. 2005 . The sustainability of western Atlantic bluefin tuna: a warm-blooded fish in a hot-blooded fishery . Bulletin of Marine Science , 76 : 363 – 384 . OpenURL Placeholder Text WorldCat Porch C. E. , Hanke A. 2018 . Estimating the fraction of western Atlantic bluefin tuna that spawn by age from size frequency data collected on the Gulf of Mexico spawning grounds . ICCAT Collective Volume of Scientific Papers , 74 : 3224 – 3233 . OpenURL Placeholder Text WorldCat Porch C. E. , Kleiber P., Turner S., Sibert J., Bailey R., Cort J. L. 1998 . The efficacy of VPA models in the presence of complicated movement patterns . ICCAT Collective Volume of Scientific Papers , 50 : 591 – 622 . OpenURL Placeholder Text WorldCat Porch C. E. , Restrepo V. R., Turner S. C., Scott G. P. 1995 . Virtual population analysis of Atlantic bluefin tuna incorporating movement and tagging data . ICCAT Collective Volume of Scientific Papers , 44 : 183 – 190 . OpenURL Placeholder Text WorldCat Porch C. E. , Turner S. C., Powers J. E. 2001 . Virtual population analysis of Atlantic bluefin tuna with alternative models of transatlantic migration: 1970–1997 . ICCAT Collective Volume of Scientific Papers , 52 : 1022 – 1045 . OpenURL Placeholder Text WorldCat Puncher G. N. , Cariani A., Maes G. E., Van Houdt J., Herten K., Cannas R., Rodriguez-Ezpeleta N., et al. 2018 . Spatial dynamics and mixing of bluefin tuna in the Atlantic Ocean and Mediterranean Sea revealed using next-generation sequencing . Molecular Ecology Resources , 18 : 620 – 638 . Google Scholar Crossref Search ADS PubMed WorldCat Punt A. E. 2019 a. Modelling recruitment in a spatial context: a review of current approaches, simulation evaluation of options, and suggestions for best practices . Fisheries Research , 217 : 140 – 155 . Google Scholar Crossref Search ADS WorldCat Punt A. E. 2019 b. Spatial stock assessment methods: a viewpoint on current issues and assumptions . Fisheries Research , 213 : 132 – 143 . Google Scholar Crossref Search ADS WorldCat Punt A. E. , Butterworth D. S., de Moor C. L., De Oliveira J. A. A., Haddon M. 2016 . Management strategy evaluation: best practices . Fish and Fisheries , 17 : 303 – 334 . Google Scholar Crossref Search ADS WorldCat Punt A. E. , Haddon M., Tuck G. N. 2015 . Which assessment configurations perform best in the face of spatial heterogeneity in fishing mortality, growth and recruitment? A case study based on pink ling in Australia . Fisheries Research , 168 : 85 – 99 . Google Scholar Crossref Search ADS WorldCat Ricker W. E. 1958 . Maximum sustained yields from fluctuating environments and mixed stocks . Journal of the Fisheries Research Board of Canada , 15 : 991 – 1006 . Google Scholar Crossref Search ADS WorldCat Rooker J. R. , Secor D. H., DeMetrio G., Kaufman A. J., Belmonte Ríos A., Tičina V. 2008 a. Evidence of trans-Atlantic movement and natal homing of bluefin tuna from stable isotopes in otoliths . Marine Ecology Progress Series , 368 : 231 – 239 . Google Scholar Crossref Search ADS WorldCat Rooker J. R. , Secor D. H., De Metrio G., Schloesser R., Block B. A., Neilson J. D. 2008 b. Natal homing and connectivity in Atlantic bluefin tuna populations . Science , 322 : 742 – 744 . Google Scholar Crossref Search ADS PubMed WorldCat Rooker J. R. , Arrizabalaga H., Fraile I., Secor D. H., Dettman D. L., Abid N., Addis P., et al. 2014 . Crossing the line: migratory and homing behaviors of Atlantic bluefin tuna . Marine Ecology Progress Series , 504 : 265 – 276 . Google Scholar Crossref Search ADS WorldCat Secor D. H. 2014 . The unit stock concept: bounded fish and fisheries. In Stock Identification Methods: Applications in Fishery Science , 2nd edn, pp. 7 – 28 . Ed. by Cadrin S. X., Kerr L. A., Mariani S.. Elsevier , Waltham . 588 pp. Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Siskey M. R. , Wilberg M. J., Allman R. J., Barnett B. K., Secor D. H. 2016 . Forty years of fishing: changes in age structure and stock mixing in northwestern Atlantic bluefin tuna (Thunnus thynnus) associated with the size-selective and long-term exploitation . ICES Journal of Marine Science , 73 : 2518 – 2528 . Google Scholar Crossref Search ADS WorldCat Szuwalski C. S. , Punt A. E. 2015 . Can an aggregate assessment reflect the dynamics of a spatially structured stock? Snow crab in the eastern Bering Sea as a case study . Fisheries Research , 164 : 135 – 142 . Google Scholar Crossref Search ADS WorldCat Taylor N. G. , McAllister M. K., Lawson G. L., Carruthers T., Block B. A. 2011 . Atlantic bluefin tuna: a novel multistock spatial model for assessing population biomass . PLoS ONE , 6 : e27693 . Google Scholar Crossref Search ADS PubMed WorldCat Thorrold S. R. , Latkoczy C., Swart P. K., Jones C. M. 2001 . Natal homing in a marine fish metapopulation . Science , 291 : 297 – 299 . Google Scholar Crossref Search ADS PubMed WorldCat Vincent M. T. , Brenden T. O., Bence J. R. 2017 . Simulation testing the robustness of a multiregion, tag-integrated assessment model that exhibits natal homing and estimates natural mortality and reporting rate . Canadian Journal of Fisheries and Aquatic Sciences , 74 : 1930 – 1949 . Google Scholar Crossref Search ADS WorldCat Webster D. G. 2009 . Adaptive Governance . The MIT Press , Cambridge . 416 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Weston A. E. , Kerr L. A., Cadrin S. X., Morse M. R. 2019 . Performance of a F0.1 management procedure using alternative operating models . ICCAT Collective Volume of Scientific Papers , 76 : 71 – 99 . OpenURL Placeholder Text WorldCat Wiedenmann J. , Wilberg M. J., Sylvia A., Miller T. J. 2015 . Autocorrelated error in stock assessment estimates: implications for management strategy evaluation . Fisheries Research , 172 : 325 – 334 . Google Scholar Crossref Search ADS WorldCat Ying Y. , Chen Y., Lin L., Gao T. 2011 . Risks of ignoring fish population spatial structure in fisheries management . Canadian Journal of Fisheries and Aquatic Sciences , 68 : 2101 – 2120 . Google Scholar Crossref Search ADS WorldCat Zarrad R. , Walter J., Lauretta M. 2018 . VPA-2BOX model diagnostics used in the 2014 assessment of eastern Atlantic bluefin tuna . ICCAT Collective Volume of Scientific Papers , 74 : 2809 – 2831 . OpenURL Placeholder Text WorldCat © International Council for the Exploration of the Sea 2020. All rights reserved. For permissions, please email: [email protected] This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Neuenfeldt,, Stefan;Bartolino,, Valerio;Orio,, Alessandro;Andersen, Ken, H;Andersen, Niels, G;Niiranen,, Susa;Bergström,, Ulf;Ustups,, Didzis;Kulatska,, Nataliia;Casini,, Michele

doi: 10.1093/icesjms/fsaa047pmid: N/A

We thank Keith Brander (KB) for his interest in, and comments (Brander, 2020) on, our paper (Neuenfeldt et al., 2020). We agree with him that increased hypoxia has played a major role in the reduced growth of Eastern Baltic cod, Gadus morhua, especially during the last two decades, and we acknowledge this in our paper. We also agree with KB that the actual mechanisms remain to be fully disclosed (Brander, 2020). As KB correctly mentions (Brander, 2020), we (Neuenfeldt et al., 2020) propose that the principal mechanism linking hypoxia and juvenile cod growth in the Eastern Baltic Sea is that increased extent of hypoxia has reduced the abundance of benthic prey suitable for post-settlement cod. We further propose that the growth reduction has resulted in a change in cod size distribution by an accumulation of smaller cod, which intensifies the competition for food in the benthivorous life stage and leads to a feedback loop that further reduces the growth. KB, on the other hand, argues that the observed changes in the oxygen environment alone may be sufficient to explain the decline in food uptake observed in our paper, because cod feed less when oxygen saturation is low (Chabot and Dutil, 1999). An implicit conclusion of KB’s arguments is then that benthic prey does not necessarily have to be scarce, as assumed in our paper, because hypoxia itself may be limiting the food intake. In this reply, we argue that the part of the benthic fauna preyed upon by cod has decreased, because this is important for our argument that small cod are food limited and hence increased density of smaller sizes due to decreased growth leads to the aforementioned feedback loop. Furthermore, we deal with the main assumption that KB bases his argument upon, namely that small cod may dwell in hypoxic waters to an extent that significantly reduces their ability to process the food. We argue that small cod probably are not subjected to a critical lack of oxygen that would decrease their physiological potential for consuming and processing food. We argue that cod feeding level has declined because of a hypoxia-related decline in benthic prey species, assuming that Saduria entomon abundance has decreased in the hypoxic parts of the seabed of the Baltic Sea and that, especially for small cod, there is no other benthic species possibly substituting S. entomon in their diet. Evidence on abundance, biomass, or distribution trends for S. entomon is sparse. Karlson et al. (2019) found for the Baltic proper that S. entomon populations in the seabed of the open sea have declined. HELCOM’s fact sheet on S. entomon states that “The loss of oxygen in the Baltic Sea bottoms has caused some declines in the population” (https://helcom.fi/media/red%20list%20species%20information%20sheet/HELCOM-Red-List-Saduria-entomon.pdf). The mass of S. entomon in the stomachs of both small cod and larger cod has decreased (Figure 1). The larger cod would have had the opportunity to reach S. entomon in hypoxic water, if it was there (Neuenfeldt et al., 2009). The mass of “other prey” (i.e. not cod, clupeids, Mysis mixta and S. entomon) has increased only in the stomachs of cod >45 cm (Fig. 2 in Neuenfeldt et al., 2020). It would be interesting to see if the increase in “other prey” at least partially was due to alternative benthos organisms, and if so, whether these potentially would be inedible for cod <45 cm. Nevertheless, we interpret the low biomass of S. entomon in the stomachs of cod of all sizes as a clear indication that this benthic isopod indeed has decreased in abundance in the distributional area of cod simultaneously with the increase in hypoxia after the mid-1990s (Fig. 7 in Neuenfeldt et al., 2020) as explained in our discussion section. Moreover, the other prey type that declined in the cod stomachs was sprat, Sprattus sprattus (Figure 1), corresponding to a decline in sprat abundance in the area of cod distribution (Casini et al. 2016), and therefore it would be logical and coherent to suppose that also S. entomon declined in cod stomachs for a decrease in this resource. Finally, the average amounts of herring (Clupea harengus) in the stomachs of cod <30 cm did not decrease, but variability between years increased (Figure 1). The amount of nekto-benthic M. mixta in the stomachs is much lower since 2008 (Figure 1). These observations would be difficult to explain, if mere exposure to hypoxia should have limited the food uptake of cod. Figure 1. Open in new tabDownload slide Stomach contents over time in Gadus morhua by average mass, given that the specific prey is found in the stomach. Note that the listed prey groups do not comprise the complete cod diet. The time series are displayed separately for 21–30 cm cod (dark closed dots) and 31–45 cm cod (light closed dots). Figure 1. Open in new tabDownload slide Stomach contents over time in Gadus morhua by average mass, given that the specific prey is found in the stomach. Note that the listed prey groups do not comprise the complete cod diet. The time series are displayed separately for 21–30 cm cod (dark closed dots) and 31–45 cm cod (light closed dots). While hypoxia generally addresses oxygen saturation below 100%, oxygen saturation well below 100% is not critical for cod (Plante et al., 1998; Chabot and Dutil, 1999) and Eastern Baltic cod >45 cm even dive down into oxygen saturations below 20% for a few hours at a time, probably in search of food (Neuenfeldt et al., 2009). Behrens et al. (2012) have shown experimentally that after such dives, there were no indications of an oxygen debt, implying that there were no indications of a negative effect on digestive capacity (Behrens et al., 2012). Hence, on a short time scale (hours), cod >45 cm can cope with very low oxygen saturation and return to waters with sufficient oxygen to evacuate their stomachs, by swimming a few metres upwards in the stratified water column, without any effects on their feeding capability. There is some indication that also smaller cod between 20 and 30 cm show this behaviour, although we agree with KB that there is no conclusive evidence. Acoustics measurements in the Bornholm Basin from March 2002 have shown that larger cod apart from the vertical feeding raids stay some metres above the bottom (Andersen et al., 2017, Fig. 6b). This applies to small cod as well, according to cod target strength groups that are not shown in Andersen et al. (2017). Measurements from otoliths further support the hypothesis that 20–30 cm cod perform vertical migrations similar to >45 cm cod (Hüssy, 2010, Fig. 3). Here, the daily fluctuating levels of opacity during the growth season imply vertical migration in a thermal gradient. Since in the Baltic Sea the thermal gradient is often overlaying with the oxygen gradient, the observed opacity pattern implies varying oxygen environment, too. Hence, both acoustics and otolith data indicate that also small cod perform vertical migration. The vertical distance small cod would have to move to get from oxygen saturation below 28% (lethal as to Plante et al., 1998) to oxygen saturation above 65% (not affecting growth according to Chabot and Dutil, 1999) is variable, but usually only few metres. For example, in observations described in Andersen et al. (2017), oxygen saturation at the bottom was 25%, while 3.5 m above bottom oxygen saturation was ≥65% and hence uncritical for food conversion and growth as experimentally deduced in Chabot and Dutil (1999). Even if the small cod would not migrate vertically, the extrapolation of results from the experiment by Chabot and Dutil (1999) is a problematic point in the arguments of KB. Except for high temperatures well above the optimum, ∼15–16°C for cod (Andersen, 2012), the limitation of maximum sustained feeding rate at high oxygen levels is due to the intrinsic physiological performance of the fish rather than the actual content of oxygen at high saturation levels (e.g. Buentello et al., 2000; Remen et al., 2016). In the Baltic Sea, cod experience temperatures within the range of 6.23 ± 2.21°C (Righton et al., 2010). Hence, the maximum sustained food consumption rate would not be higher in the Baltic Sea because oxygen content at saturation in the colder and less saline water is higher as compared to the conditions of 10°C and salinity of 28 in the experiments of Chabot and Dutil (1999). Therefore, the relationship between oxygen saturation and maximum feeding rate estimated by Chabot and Dutil (1999) cannot be extrapolated to 100% oxygen saturation in the Baltic Sea for the provision of 40 g/d maximum feeding rate in the Baltic Sea as compared to 32.6 g/d by Chabot and Dutil (for cod with an initial size of 44.1 ± 3.1 cm total length; 1999). On the contrary, the lower temperatures in the Baltic Sea imply lower physiological rates including metabolic expenses and maximum feeding rates. The decrease in feeding level to 68% of the maximum sustained feeding rate due to low O2 concentrations of 4.73 ml/l as suggested by KB is therefore not realistic and probably heavily overestimated. We would also like to clarify that, in Limburg and Casini (2019), the analyses of the effects of low-oxygen exposure (indexed by Mn/Mg in otoliths) on cod condition were all in relative terms and not related to a specific threshold of 2 ml/l. We agree with KB that, given the problems the small cod apparently are facing currently, their behaviour and feeding biology should be a topic of high priority. Further dietary analyses based on the abundantly available cod stomach data, as well as time series of the benthos that is consumed by small cod, would be desirable. We also agree that it would be valuable to have better information on the behaviour of small cod. Further empirical knowledge on the direct effects of hypoxia on feeding is needed from experimental work, as hypoxia may be a factor contributing to lowering the ingestion rates, which in the future can be expected to be more severe as a consequence of increasing temperatures. We furthermore think that targeted field observations using hydroacoustics and data storage tags are needed to answer the questions on individual behaviour, which form the basis to disentangling the effects of hypoxia and benthic food availability on the poor condition and slow growth of Eastern Baltic cod. Neuenfeldt, S., Bartolino, V., Orio, A., Andersen, K. H., Andersen, N. G., Niiranen, S., Bergström, U., Ustups, D., Kulatska, N., and Casini, M. 2020. Reply to “Reduced growth in Baltic Sea cod may be due to mild hypoxia”—a comment to Neuenfeldt et al. (2020). – ICES Journal of Marine Science, 00:000–000. References Andersen N. G. 2012 . Influences of potential predictor variables on gastric evacuation in Atlantic cod Gadus morhua feeding on fish prey: parameterization of a generic model . Journal of Fish Biology , 80 : 595 – 612 . Google Scholar Crossref Search ADS PubMed WorldCat Andersen N. G. , Lundgren B., Neuenfeldt S., Beyer J. E. 2017 . Diel vertical interactions between Atlantic cod Gadus morhua and sprat Sprattus sprattus in a stratified water column . Marine Ecology Progress Series , 583 : 195 – 209 . Google Scholar Crossref Search ADS WorldCat Behrens J. W. , Axelsson M., Neuenfeldt S., Seth H. 2012 . Effects of hypoxic exposure during feeding on SDA and postprandial cardiovascular physiology in the Atlantic cod, Gadus morhua . PLoS One , 7 : e46227 . Google Scholar Crossref Search ADS PubMed WorldCat Brander K. 2020 . Reduced growth in Baltic Sea cod may be due to mild Hypoxia—a comment to Neuenfeldt et al. (2019) . ICES Journal of Marine Science, doi:10.1093/icesjms/fsaa041 . OpenURL Placeholder Text WorldCat Buentello J. A. , Gatlin D. M.III Neill W. H. 2000 . Effects of water temperature and dissolved oxygen on daily feed consumption, feed utilization and growth of channel catfish (Ictalurus punctatus) . Aquaculture , 182 : 339 – 352 . Google Scholar Crossref Search ADS WorldCat Casini M. , Käll F., Hansson M., Plikshs M., Baranova T., Karlsson O., Lundström K. et al. 2016 . Hypoxic areas, density-dependence and food limitation drive the body condition of a heavily exploited marine fish predator . Royal Society Open Science , 3 : 160416 . Google Scholar Crossref Search ADS PubMed WorldCat Chabot D. , Dutil J. D. 1999 . Reduced growth of Atlantic cod in non-lethal hypoxic conditions . Journal of Fish Biology , 55 : 472 – 491 . Google Scholar Crossref Search ADS WorldCat Hüssy K. 2010 . Why is age determination of Baltic cod (Gadus morhua) so difficult? ICES Journal of Marine Science , 67 : 1198 – 1205 . Google Scholar Crossref Search ADS WorldCat Karlson A. M. , Gorokhova E., Gårdmark A., Pekcan-Hekim Z., Casini M., Albertsson J., Sundelin B. et al. 2019 . Linking consumer physiological status to food-web structure and prey food value in the Baltic Sea . Ambio , 49 : 1 – 16 . Google Scholar PubMed OpenURL Placeholder Text WorldCat Limburg K. E. , Casini M. 2019 . Otolith chemistry indicates recent worsened Baltic cod condition is linked to hypoxia exposure . Biology Letters , 15 : 20190352 . Google Scholar Crossref Search ADS PubMed WorldCat Neuenfeldt S. , Andersen K. H., Hinrichsen H. H. 2009 . Some Atlantic cod Gadus morhua in the Baltic Sea visit hypoxic water briefly but often . Journal of Fish Biology , 75 : 290 – 294 . Google Scholar Crossref Search ADS PubMed WorldCat Neuenfeldt S. , Bartolino V., Orio A., Andersen K. H., Andersen N. G., Niiranen S., Bergström U. et al. 2020 . Feeding and growth of Atlantic cod (Gadus morhua L.) in the Eastern Baltic Sea under environmental change . ICES Journal of Marine Science , 77 : 624 – 632 . Google Scholar Crossref Search ADS WorldCat Plante S. , Chabot D., Dutil J. D. 1998 . Hypoxia tolerance in Atlantic cod . Journal of Fish Biology , 53 : 1342 – 1356 . Google Scholar Crossref Search ADS WorldCat Remen M. , Sievers M., Torgersen T., Oppedal F. 2016 . The oxygen threshold for maximal feed 172 intake of Atlantic salmon post-smolts is highly temperature-dependent . Aquaculture , 464 : 582 – 592 . Google Scholar Crossref Search ADS WorldCat Righton D. A. , Andersen K. H., Neat F., Thorsteinsson V., Steingrund P., Svedäng H., Michalsen K. et al. 2010 . Thermal niche of Atlantic cod Gadus morhua: limits, tolerance and optima . Marine Ecology Progress Series , 420 : 1 – 13 . Google Scholar Crossref Search ADS WorldCat © International Council for the Exploration of the Sea 2020. All rights reserved. For permissions, please email: [email protected] This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Denechaud,, Côme;Smoliński,, Szymon;Geffen, Audrey, J;Godiksen, Jane, A

doi: 10.1093/icesjms/fsz259pmid: N/A

Abstract Otolith shape analysis provides a robust tool for the discrimination of many fish stocks in the context of fisheries management. However, there has been little research to examine within-stock temporal stability of otolith morphology in relation to changes in the environment and stock conditions. This study investigated the variability of Northeast Arctic (NEA) cod (Gadus morhua) otolith shape between 1933 and 2015, using elliptical Fourier descriptors extracted from archived material of 2968 mature fish. Series of hierarchical multivariate models were developed to relate shape to the identified optimal windows of some environmental drivers. Differences between years accounted for <3% of the observed variation and no significant differences were found between the average cohort shapes. The models not only confirmed that fish growth was the strongest driver of shape differences, but also highlighted effects of temperature and biomass-related variables at different life stages. Extrinsic factors described only a small fraction of the observed variance, which indicates that environmental changes over time likely account for less than the natural inter-individual variability in otolith shape. These results suggest that overall shape remains relatively stable through time within NEA cod, which further contributes towards a consensus on the biological interpretation of shape differences. Introduction Otoliths are widely used in fisheries science (Campana, 2005). They form annual and daily growth structures, are metabolically inert, and hold permanent records of life-history events, making them reliable indicators of individual fish age and population age structure (Campana, 1999, 2001). Otolith shape is also specific to species and often populations within species, and shape analysis can support the identification and discrimination of fish stocks. This analysis provides a basis for understanding population dynamics, which is of primary importance for the conservation and successful management of marine fish resources (Hammer and Zimmermann, 2005). Analysis of otolith shape has been successfully applied to stock discrimination in multiple instances (Campana and Casselman, 1993; Stransky and MacLellan, 2005; Petursdottir et al., 2006; Stransky et al., 2008). While otolith traits such as annuli spacing are well understood and biologically interpretable, the determinants of otolith shape remain less clear. Inter-population variations in otolith shape have been linked to both genetic and environmental influences (Campana and Casselman, 1993; Lombarte and Lleonart, 1993; Cardinale et al., 2004; Stransky and MacLellan, 2005; Vignon and Morat, 2010), but these studies also highlighted a strong variability related to individual-specific factors such as year-class, sex or age, as well as to growth and the local environment. While overall shape is genetically constrained and species- or stock-specific (Lombarte and Lleonart, 1993), a non-negligible proportion of the observed variations seems to be of environmental origin but few analyses quantify or distinguish the relative influence of each factor (Campana and Casselman, 1993; Cardinale et al., 2004; Vignon and Morat, 2010). Local intra-population variability in otolith shape has been little studied in comparison with larger-scale variations (Vignon, 2015). While the latter is expected to be associated with significant shape differences, local environment has also been identified to be an important contributor to shape variability (Vignon, 2017). Yet, studies aiming at discriminating stocks based on otolith shape often neglect potential sources of within-stock variability in favour of the variations found among different stocks. Castonguay et al. (1991) emphasized the need to separate stocks with caution when relying on otolith morphology because year-class effects between samples could be mistaken for stock difference, and most stock discrimination studies successfully apply these recommendations. It was, therefore, suggested that the characteristics of otolith shape used in stock discrimination should be recalculated each year (Begg and Brown, 2000). However, the underlying assumption that otolith shape can vary between years has seldom been tested. Jónsdóttir et al. (2006) reported greater differences between locations than among years within a location for Icelandic cod, and Vignon (2015) reported similar findings for eight species of tropical fish from Pacific ocean islands. However, both studies used only otoliths from two or three consecutive years and to date, there has not been a comprehensive attempt at evaluating its stability through longer periods of time. The present study quantifies otolith shape variability within a single stock unit and discusses its variability in relation to environmental changes. This study investigated the temporal stability of otolith shape in Northeast arctic (NEA) cod, the stock designation for the large migratory Atlantic cod (Gadus morhua) population inhabiting the Barents Sea. NEA cod otoliths have been collected during fisheries surveys for nearly a century, and this archival collection is ideal for assessing the within-stock temporal variability of sagittal otolith shape, to test the hypothesis that otolith morphology has changed throughout the last century in response to changing environmental conditions. Otoliths from a single stock were collected from archived material and their outline described using elliptical Fourier descriptors. Series of hierarchical multivariate models were then developed to compare these descriptors between years and to relate them with different environmental variables. Conditions during a fish’s early life stages could be the most influential on otolith shape, and to verify this hypothesis a sliding window analysis was conducted to identify the optimal time windows of the predictors driving shape differences. Quantifying the stability of otolith shape over the last century could bring new insights about its resilience to changing environmental conditions, which can be of importance to understanding and identifying the factors driving shape variability. Material and methods Sample collection At the Institute of Marine Research in Bergen (Norway), otoliths have been routinely used for age determination purposes since 1932, and, for most of these, information on catch location, date, and biological parameters (age, length, weight, and sex) are available. Sagittal otoliths of mature NEA cod from the period 1933 to 2015 were randomly selected and retrieved from the extensive archive available at the institute. The otoliths had all been previously broken for age estimation, and so only one whole otolith from each pair could be collected intact. Consequently, the samples studied here consisted of a mixture of left and right otoliths from different fish, without the possibility to include both from the same pair. Norwegian coastal cod (NCC) is the stock designation used to refer to a collection of genetically distinct local coastal populations that exhibit different life-history traits and environmental exposures from the migratory NEA cod population (Salvanes et al., 2004). Both stocks have been successfully separated using inner zone morphology, otolith outer shape and genetic methods (Rollefsen, 1934; Berg et al., 2005; Stransky et al., 2008; Dahle et al., 2018; Johansen et al., 2018). To prevent the inclusion of otoliths from NCC in the present study, only fish classified “certain NEA cod” by experienced age readers based on inner otolith morphology (Stransky et al., 2008) were selected. The mean annual proportion of uncertain NEA cod otolith in the archives was only 5.24%. Proportions of otoliths assigned to the different cod stock in the Barents Sea ecoregion are presented in detail in Supplementary material 1. All the samples used in this analysis were fished in the Lofoten spawning grounds (Figure 1), which encompasses a large area along the north-western Norwegian coast, comprising both the Lofoten archipelago and the Vesterålen region. Because there was no single source going back as far as 1933, this selection comprised a mixture of survey-caught fish as well as samples collected from commercial landings. Preliminary analysis showed a potential selectivity bias of gillnets towards bigger, faster-growing fish, and samples caught that way were thus excluded from the selection, limiting it to either bottom trawls, longlines, or seine. Figure 1. Open in new tabDownload slide Map of the Northeast Atlantic region, including the northern coast of Norway and the Barents Sea with added depth contour. Framed area represents the Lofoten spawning grounds. Blue line represents the Kola transect of in situ temperature measurements. Red arrows represent the general direction of NEA cod spawning migration. Figure 1. Open in new tabDownload slide Map of the Northeast Atlantic region, including the northern coast of Norway and the Barents Sea with added depth contour. Framed area represents the Lofoten spawning grounds. Blue line represents the Kola transect of in situ temperature measurements. Red arrows represent the general direction of NEA cod spawning migration. The initial sampling comprised >4000 individuals of age 7–21 of which 3913 otoliths were intact. For the purposes of this study, only fish of age 8 were retained to limit the influence of age on the shape analysis and focus on between-cohort (year-class) differences. This yielded a total of 2968 otoliths and at least 10 otoliths per year, representing 1467 females, 1388 males, and 113 fish of unspecified sex, divided between 1456 left-hand and 1512 right-hand orientations (see Supplementary material 1 for detailed information about the repartition of samples per year). Otolith shape analysis Each otolith was first manually cleaned to remove residual tissue and then weighed on a high-resolution scale. Individual images of the whole otolith were taken under reflected light using a Nikon SMZ25 stereomicroscope mounted with a Nikon Digital Sight DS-Fi2 camera and a Nikon SHR Plan Apo 0.5× WD:71 lens. The otoliths were all positioned with the proximal surface and sulcus acusticus facing up, and the dorsal side at the top of the image. Acquisition parameters such as shutter speed, aperture, white balance and sensitivity (ISO) were kept uniform between all samples. Images were processed for outline extraction directly in R (R Core Team, 2016) with the packages imager (Barthelme, 2019) and Momocs (Bonhomme et al., 2014). Each image was first converted to greyscale, then colour-inverted and binarized to generate a dark shape on a white background. Because the otoliths comprised both orientations, right otoliths were mirrored horizontally during this processing step. The outlines were then detected and extracted by intensity thresholding, based on the transition from black to white pixel values. Because of fibre leftovers or black mould spots caused by humidity, the outside image boundary of some of the older otoliths was often deteriorated. The first batch of outlines generated was consequently superimposed on their associated original pictures and each of them was visually inspected for failures or artefacts in the detected outline. Whenever the superimposed outline deviated from the otolith shape, the original image was imported into Adobe Photoshop 2019, corrected, and reprocessed. To analyse the variations in shape, an elliptical Fourier analysis (EFA) was performed on the delineated otolith contours (Lestrel, 1997). While cod otoliths often display finer scale crenulations and lobes that require high numbers of harmonics to be recreated with precision, too many will often come at the cost of computational speed. For each otolith, the first 99 harmonics (H) were thus arbitrarily extracted and normalized by the first to provide elliptical Fourier descriptors (EFDs) invariant to size, rotation, or starting point (Kuhl and Giardina, 1982). In addition, the extraction path was made homologous for all otoliths by selection of a starting point located to the left of the centroid. This was done to further prevent mirroring effects wherein unguided outline extraction can sometimes produce inverted shapes (Bonhomme et al., 2014). To determine the number of harmonics needed to optimally reconstruct the otolith contour, the cumulated Fourier power (PF; Lestrel, 1997) was calculated on a subset of 100 otoliths randomly selected across the whole dataset. The number of harmonics used in the analysis (nmax) was then chosen such that PF (nmax) explained 99.99% of the variance in contour coordinates, i.e. 99.99% of the otolith shape obtained at 99 harmonics was recreated at nmax. Statistical analyses To investigate the relation between shape and environment, a shape matrix S was built with individual fish in rows and corresponding EFDs in columns. A total of three explanatory matrices were evaluated in relation to the shape matrix S: A year matrix Y, to investigate potential variations of shape over time. Since every fish was of the same age, the years were used as proxies for cohorts (or year class) and likely reflected inter-annual differences. An individual matrix I, grouping biological variables related to each sample. This matrix was included as conditional factor to remove the effects of individual variables on otolith shape. It was composed of fish length (L) to account for shape variation related to growth rates differences between individuals, sex (Sx) as a factor potentially affecting fish metabolism and otolith mineralization, and otolith orientation (Or) to account for potential cofounding effects of directional asymmetry (Mahé et al., 2019). An environmental matrix E, grouping external variables related to the environment and population dynamics, composed of four main variables. The monthly and yearly Barents Sea temperatures (T°) of the 0–200 m depth layer between 1921 and 2015 were extracted from the Kola section in situ records provided by the Russian Polar Research Institute of Marine Fisheries and Oceanography (Bochkov, 1982; Tereshchenko, 1996). This proxy offers a good representation of climatological conditions within the area occupied by the NEA cod both during its juvenile and adult stages (Dippner and Ottersen, 2001). The total stock biomass (TS), spawning stock biomass (SSB), and recruitment (R) were obtained from the Arctic Fisheries Working Group at the International Council for the Exploration of the Sea for the period 1946–2015 (ICES, 2018), and from the extended virtual population analysis carried out by Hylen (2002) for the period 1920–1945. Hierarchical redundancy analyses (RDA) were then conducted to model the matrix S in relation to the different explanatory matrices with a statistical protocol similar to Mille et al. (2016). RDA is an extension of principal component analysis, wherein the variation in a set of response variables is related to a set of explanatory variables through the extension of multiple linear regressions to multivariate data (Legendre and Legendre, 2012). For each model, the exploratory matrix of interest was first tested in relation to S through a normal RDA. Then, a similar partial RDA (pRDA) was built using the matrix I as a conditional factor, which removed the variance in shape associated with fish-specific variables before constraining with the tested explanatory matrix. This step was done to ensure that potential significant differences through time or environmental effects found in the initial RDA were not caused by any measured confounding factors, in case of a strong joint effect. To test the significance of each model or explanatory matrix, ANOVA permutation tests (999 iterations, marginal effects, Type II) were performed (Legendre and Legendre, 2012). In addition, the potential collinearity between explanatory variables was investigated by calculating their variance inflation factors (VIF) with a more conservative threshold defined at VIF <2 (Borcard et al., 2018). No evidence for multicollinearity was found in any of the analyses. Finally, variation partitioning was performed to quantify the respective contributions of Y, E, and I to otolith shape differences within each model. The environmental descriptors associated with each cohort usually describe only a snapshot of an older individual’s exposure because they lack fine temporal resolution, especially for pre-instrument records. Such environmental proxies are not representative of the diversity of environmental conditions experienced throughout the life of an organism, nor will they reveal the contrasted effects similar changes can have at different times of the year (Cook et al., 2012; Kruuk et al., 2015; Roberts et al., 2015). This is especially true for mobile organisms like NEA cod, which starts migrating thousands of kilometres southward to the spawning grounds once it reaches sexual maturity, and thus experiences diverse conditions throughout its life cycle (Robichaud and Rose, 2004). A preliminary analysis sequentially compared models built using values of environmental variables at the year of catch to those at the year of hatch and found significant differences between their respective contribution to shape differences (see Supplementary material 2). To account for this problem, optimal time windows of the four environmental variables were explored with a sliding window analysis, which allowed for the statistical identification of the best predictors and their associated critical time windows (van de Pol et al., 2016). Different time windows were investigated sequentially within an 11-year period (counted backwards from the year of capture), covering both the life of the fish and the conditions up to 2 years prior to its hatch in case of eventual lagged effects. Collinearity of identified environmental predictors was tested with Pearson’s correlation tests. Because the high number of time windows analysed could potentially lead to the misidentification of an environmental signal (or false positive), Monte Carlo randomization tests were then conducted to assess the probability of obtaining similar strong statistical support of the model by chance (van de Pol et al., 2016; Smoliński, 2019). Detailed statistical methodology is available in Supplementary material 3. The illustrated workflow is presented in Figure 2. Figure 2. Open in new tabDownload slide Schematic representation of the hierarchical RDA analyses performed in the EFA: (1) the global relationship between shape and cohort (year-class), (2) the influence of environmental conditions at the identified optimal windows, and (3) the global influence of both environmental and biological variables. Steps (a) and (b), respectively, refer to the RDA and associated partial RDA (where the variance attributed to matrix I was removed). Figure 2. Open in new tabDownload slide Schematic representation of the hierarchical RDA analyses performed in the EFA: (1) the global relationship between shape and cohort (year-class), (2) the influence of environmental conditions at the identified optimal windows, and (3) the global influence of both environmental and biological variables. Steps (a) and (b), respectively, refer to the RDA and associated partial RDA (where the variance attributed to matrix I was removed). This first analysis utilized a large number of individual shapes, but because the environmental values associated with each fish were consistently repeated within the same cohort, the identification of an environmental signal was potentially masked by inter-individual variability. In a second analysis, the sliding window analysis and environmental models were thus recreated using a shape matrix Syear composed of EFDs averaged per year, wherein each cohort (N = 78) was represented by a single average shape. Similarly, the individual matrix I was also replaced by Iyear, where the mean fish length for each year was computed as the mean cohort body length at age 8. This cohort-based analysis achieved three aims: to reduce variability associated with between-individual shape differences, to determine if average shape changed during the period investigated, and to directly relate it to environmental factors. All analyses were conducted using the R scientific computing language (R Core Team, 2016) and following packages: vegan (Oksanen et al., 2019), PerformanceAnalytics (Peterson et al., 2019), and tidyverse (Wickham, 2017). Results Analysis of individual shapes Based on the Fourier power calculated at 99 harmonics, the first 54 harmonics explained at least 99.99% of the variation in otolith contour (Figure 3) and were thus selected for further analysis. Year-to-year variations in shape were significant but only accounted for a minor portion of the observed variability. In models 1.a and 1.b (respectively, Y only and Y with conditional removal of I), matrix Y explained 2.67% and 2.41% of the variability in otolith shape (Table 1). Variation partitioning (Figure 4) indicated that the matrix I used as a set of conditional factors explained around 2.2% of the variability in shape, with a joint contribution of both matrices equal to 0.3%. Figure 3. Open in new tabDownload slide Fourier reconstruction of an NEA cod otolith at 1 (a), 5 (b), and 54 (c) harmonics. Corresponding number of elliptical Fourier descriptors after normalization of the first harmonic is indicated. Figure 3. Open in new tabDownload slide Fourier reconstruction of an NEA cod otolith at 1 (a), 5 (b), and 54 (c) harmonics. Corresponding number of elliptical Fourier descriptors after normalization of the first harmonic is indicated. Figure 4. Open in new tabDownload slide Venn diagrams of variation partitioning between the matrices (year Y, individual I, and environmental Eopti) tested in the analysis of individual shapes (a, b) and the analysis of average cohort shapes (c). Values in the non-overlapping parts of each circle represent the strict contribution of the corresponding matrix to the model. The value in the overlapping section represents the joint contribution of both matrices. Figure 4. Open in new tabDownload slide Venn diagrams of variation partitioning between the matrices (year Y, individual I, and environmental Eopti) tested in the analysis of individual shapes (a, b) and the analysis of average cohort shapes (c). Values in the non-overlapping parts of each circle represent the strict contribution of the corresponding matrix to the model. The value in the overlapping section represents the joint contribution of both matrices. Table 1. Results of the hierarchical RDA models for individual shapes detailed in Figure 2. . Variables . Model . Year . T° . TS . SSB . L . Sx . Or . df . % . p-Value . Model 1 Year (a) 5.06% 0.001 73 2.67 0.001 (b) 4.86% 0.001 Conditional 73 2.41 0.001 Model 2 Optimal windows (a) 0.14% 0.001 0.80% 0.001 0.42% 0.001 3 1.15% 0.001 (b) 0.14% 0.003 0.89% 0.001 0.32% 0.001 Conditional 3 1.11 0.001 Model 3 Global 0.14% 0.001 0.87% 0.001 0.31% 0.001 1.33% 0.001 0.15% 0.010 0.94% 0.001 7 3.54 0.001 . Variables . Model . Year . T° . TS . SSB . L . Sx . Or . df . % . p-Value . Model 1 Year (a) 5.06% 0.001 73 2.67 0.001 (b) 4.86% 0.001 Conditional 73 2.41 0.001 Model 2 Optimal windows (a) 0.14% 0.001 0.80% 0.001 0.42% 0.001 3 1.15% 0.001 (b) 0.14% 0.003 0.89% 0.001 0.32% 0.001 Conditional 3 1.11 0.001 Model 3 Global 0.14% 0.001 0.87% 0.001 0.31% 0.001 1.33% 0.001 0.15% 0.010 0.94% 0.001 7 3.54 0.001 Adjusted R2 for partial RDA is calculated using Ezekiel’s formula as the fraction of variance explained by that model after removing the fraction associated with the eventual covariables. “Variables” gives the percentage of explained variance and p-values for each significant variable within a model (p < 0.05). “Model” gives the df and the adjusted percentage of variation (based on adjusted R2) explained by each model with associated p-value. df, degrees of freedom. Open in new tab Table 1. Results of the hierarchical RDA models for individual shapes detailed in Figure 2. . Variables . Model . Year . T° . TS . SSB . L . Sx . Or . df . % . p-Value . Model 1 Year (a) 5.06% 0.001 73 2.67 0.001 (b) 4.86% 0.001 Conditional 73 2.41 0.001 Model 2 Optimal windows (a) 0.14% 0.001 0.80% 0.001 0.42% 0.001 3 1.15% 0.001 (b) 0.14% 0.003 0.89% 0.001 0.32% 0.001 Conditional 3 1.11 0.001 Model 3 Global 0.14% 0.001 0.87% 0.001 0.31% 0.001 1.33% 0.001 0.15% 0.010 0.94% 0.001 7 3.54 0.001 . Variables . Model . Year . T° . TS . SSB . L . Sx . Or . df . % . p-Value . Model 1 Year (a) 5.06% 0.001 73 2.67 0.001 (b) 4.86% 0.001 Conditional 73 2.41 0.001 Model 2 Optimal windows (a) 0.14% 0.001 0.80% 0.001 0.42% 0.001 3 1.15% 0.001 (b) 0.14% 0.003 0.89% 0.001 0.32% 0.001 Conditional 3 1.11 0.001 Model 3 Global 0.14% 0.001 0.87% 0.001 0.31% 0.001 1.33% 0.001 0.15% 0.010 0.94% 0.001 7 3.54 0.001 Adjusted R2 for partial RDA is calculated using Ezekiel’s formula as the fraction of variance explained by that model after removing the fraction associated with the eventual covariables. “Variables” gives the percentage of explained variance and p-values for each significant variable within a model (p < 0.05). “Model” gives the df and the adjusted percentage of variation (based on adjusted R2) explained by each model with associated p-value. df, degrees of freedom. Open in new tab The first iteration of the sliding window analysis for environmental variables identified total stock from the second to seventh year of fish life as the optimal environmental signal (Figure 5). Similarly, the next iteration identified mean SSB from the third to seventh year of fish life as the second-best predictor for shape variation between the different variables. Both predictors showed a wide critical window covering late juvenile stages to average age of sexual maturity. During the third and final iteration, mean temperature from September to December of the year prior to hatch was identified as the third best predictor. Contrary to the first two iterations, this step revealed a narrower period strongly supported by models, although the overall window was generally situated from hatch to recruitment age. Monte Carlo tests conducted with 1000 iterations (Figure 5) indicated that the likelihood of obtaining similar strong signals by chance was minimal for total stock and SSB (p < 0.001). The probability obtained for the temperature signal was higher but still significantly different to chance (p = 0.045). Correlation tests conducted on the identified optimal predictors revealed that total stock and SSB had a moderate positive correlation (R = 0.47, p < 0.001). The temperature had a low positive correlation with SSB (R = 0.29, p < 0.001) but no correlation with total stock. Because recruitment was not identified as an optimal predictor during the sliding window analysis and had no significant effect in any of the preliminary models, the variable was removed from matrix E. Figure 5. Open in new tabDownload slide Results of the sliding window analysis on individual shapes for: (a) total stock (first best predictor), (b) SSB (second-best predictor), and (c) temperature (third best predictor). Squares represent time windows with corresponding years (a, b) and months (c) of fish life at which it was opened or closed. Axes are graduated with years for readability. Intersecting bold lines correspond to theoretical hatch of an individual. Colour gradient shows the deviance of model generated for each window where lowest deviances indicate best model fit. The green circle and dotted lines indicate the identified optimal window. Corresponding histograms of randomized deviance distributions (1000 iterations) are shown in (2) for: (d) total stock, (e) SSB, and (f) temperature. Dashed lines indicate the deviance of the optimal window and p-values the probability of obtaining the same signal by chance. Figure 5. Open in new tabDownload slide Results of the sliding window analysis on individual shapes for: (a) total stock (first best predictor), (b) SSB (second-best predictor), and (c) temperature (third best predictor). Squares represent time windows with corresponding years (a, b) and months (c) of fish life at which it was opened or closed. Axes are graduated with years for readability. Intersecting bold lines correspond to theoretical hatch of an individual. Colour gradient shows the deviance of model generated for each window where lowest deviances indicate best model fit. The green circle and dotted lines indicate the identified optimal window. Corresponding histograms of randomized deviance distributions (1000 iterations) are shown in (2) for: (d) total stock, (e) SSB, and (f) temperature. Dashed lines indicate the deviance of the optimal window and p-values the probability of obtaining the same signal by chance. Using the identified optimal environmental windows, models 2.a and 2.b explained 1.15% and 1.11% of shape variability (Table 1). Total stock and SSB were the most significant variables (p < 0.001) and explained the highest percentage of variation (respectively, 0.80%/0.42% in model a and 0.89%/0.32% in model b). The temperature was slightly less significant in model 2.b (p = 0.003) and explained only 0.14% of shape variation in both models. When matrix I was included as an exploratory matrix, model 3 explained 3.54% of the total variation in otolith shape. Both body length and otolith orientation had the most significant effect on shape (p < 0.001) and explained, respectively, 1.33% and 0.94% of the variability. Fish sex was less significant (p = 0.010) and only accounted for 0.15%. Matrix E had a similar contribution as in model 2.b. The variation partitioning (Figure 4) revealed no joint contribution of matrices Eopti and I. Analysis of average cohort shapes Visual exploration of the morphospace occupation prior to analysis showed no significant clustering between individuals from different cohorts (Figure 6). When visualized, average cohort shapes reconstructed with the inverse Fourier transformation showed a consistent overlap (Figure 6), further demonstrating the low contribution of year-to-year variations previously found and the overall temporal stability of otolith shape. Figure 6. Open in new tabDownload slide Morphospace occupation of individual (a) and average cohort (b) otolith contours. Each cohort is associated with a colour. Figure 6. Open in new tabDownload slide Morphospace occupation of individual (a) and average cohort (b) otolith contours. Each cohort is associated with a colour. The previous statistical analysis was then repeated using these average cohort shapes. During the second sliding window analysis, total stock was again identified as the first optimal best predictor in the first iteration, with a critical window from third to sixth year of fish life. On the second iteration, the identified second-best predictor was SSB from the third to seventh year of fish life, which was identical to the second-best predictor identified in during the sliding window analysis for individual shapes. The final iteration identified mean temperature during the end of second year of fish life as the third optimal predictor (Supplementary material 1, Figure S4). Correlation tests carried out on the identified predictors revealed moderate correlation between total stock and SSB (R = 0.49, p < 0.001). However, correlation between temperature and both total stock and SSB was moderate but higher than in the first sliding window analysis (R = 0.27, p = 0.01 and R = 0.38, p < 0.001). Monte Carlo tests conducted at 1000 iterations found again a minimal probability of finding these optimal windows for total stock and SSB by chance (p = 0.002 and p = 0.005), confirming the strength of the signal. However, the probability of the identified critical windows for temperature to be a product of chance was significantly higher (p = 1.505) and above the generally accepted threshold (p = 0.05). The predictor was still retained and included in the following redundancy analysis to see if they contributed significantly to the models. Using the identified optimal environmental windows, models 2.a and 2.b explained 20.03% and 17.00% of shape variability, respectively (Table 2). The total stock was the main significant variable in both models (p < 0.001) and explained, respectively, 13.99% and 15.63% of shape variability. The contribution of SSB was also significant but lower (8.59%, p = 0.001 and 5.75%, p = 0.001). While the identified optimal window for temperature showed a higher probability to be a product of chance, the variable still contributed significantly to the models (respectively, 2.61% and 2.73%) but was less significant (p = 0.037 and p = 0.031). When matrix I was included as an exploratory matrix, model 3 explained 26.77% of the total variation in otolith shape. Average cohort body length had a significant effect on shape (p < 0.001) and explained 7.36% of its variability. Matrix Eopti had a similar contribution as in model 2.b. The variation partitioning (Figure 4) revealed a significant joint contribution of Eopti and I (3.3%), indicating that some of the environmental effects on otolith shape take place indirectly through changes in fish growth. Table 2. Results of the hierarchical RDA models for average cohort shapes. . Variables . Model . T° . TS . SSB . L . df . % . p-value . Model 2 Optimal windows (a) 2.61% 0.037 13.99% 0.001 8.59% 0.001 3 20.03 0.001 (b) 2.73% 0.031 15.63% 0.001 5.75% 0.001 Conditional 3 17.00 0.001 Model 3 Global 2.53% 0.033 14.48% 0.001 5.33% 0.002 7.36% 0.001 4 26.77 0.001 . Variables . Model . T° . TS . SSB . L . df . % . p-value . Model 2 Optimal windows (a) 2.61% 0.037 13.99% 0.001 8.59% 0.001 3 20.03 0.001 (b) 2.73% 0.031 15.63% 0.001 5.75% 0.001 Conditional 3 17.00 0.001 Model 3 Global 2.53% 0.033 14.48% 0.001 5.33% 0.002 7.36% 0.001 4 26.77 0.001 Adjusted R2 for partial RDA is calculated using Ezekiel’s formula as the fraction of variance explained by that model after removing the fraction associated with the eventual covariables. “Variables” gives the percentage of explained variance and p-values for each significant variable within a model (p < 0.05, dash indicates non-significance). “Model” gives the df and the adjusted percentage of variation (based on adjusted R2) explained by each model with associated p-values. df, degrees of freedom. Open in new tab Table 2. Results of the hierarchical RDA models for average cohort shapes. . Variables . Model . T° . TS . SSB . L . df . % . p-value . Model 2 Optimal windows (a) 2.61% 0.037 13.99% 0.001 8.59% 0.001 3 20.03 0.001 (b) 2.73% 0.031 15.63% 0.001 5.75% 0.001 Conditional 3 17.00 0.001 Model 3 Global 2.53% 0.033 14.48% 0.001 5.33% 0.002 7.36% 0.001 4 26.77 0.001 . Variables . Model . T° . TS . SSB . L . df . % . p-value . Model 2 Optimal windows (a) 2.61% 0.037 13.99% 0.001 8.59% 0.001 3 20.03 0.001 (b) 2.73% 0.031 15.63% 0.001 5.75% 0.001 Conditional 3 17.00 0.001 Model 3 Global 2.53% 0.033 14.48% 0.001 5.33% 0.002 7.36% 0.001 4 26.77 0.001 Adjusted R2 for partial RDA is calculated using Ezekiel’s formula as the fraction of variance explained by that model after removing the fraction associated with the eventual covariables. “Variables” gives the percentage of explained variance and p-values for each significant variable within a model (p < 0.05, dash indicates non-significance). “Model” gives the df and the adjusted percentage of variation (based on adjusted R2) explained by each model with associated p-values. df, degrees of freedom. Open in new tab Discussion Stability of NEA cod otolith morphology Earlier studies have acknowledged the risks of mixing individuals from different year-classes and ages when using otolith shape as a stock discrimination tool. Begg and Brown (2000) warned that a discriminant analysis of shape variables could incorrectly separate a population based on significant differences, which were related in fact to inter-individual variability. Although they suggested that the baselines for stock separation be reconsidered each year, otolith morphology studies have often assumed that the variations of shape within a particular stock were minimal through time. The present study found that the year-class effect, although significant, accounted for <3% of the overall variation in NEA cod otolith shape over an 80-year period. Because all individuals included in the analysis were of the same age-class (8 years old), this study also minimized differences related to age or maturity status effects, previously identified as strong sources of variability (Cardinale et al., 2004; Hüssy, 2008; Capoccioni et al., 2011). The weak effect found thus corroborates the common assumption that shape remains largely consistent through time within a single population. Within-stock shape stability should now be studied for a broader range of species to determine whether temporal stability can always be assumed within a stock, or if stability itself is species-specific. Studying the temporal changes of shape within different stocks from the same species, where genetic variations are limited, would also provide essential insight on how much fish somatic and otolith growth is influenced by environmental parameters and under what conditions. Some discrimination studies based on otolith shape have used individuals from multiple age- and/or year-classes (Campana and Casselman, 1993; Friedland and Reddin, 1994; Jónsdóttir et al., 2006; Tracey et al., 2006; Mahe et al., 2016) and found no significant differences between years of sampling. While this study focused on evaluating the temporal stability of otolith shape within a single stock, it is also essential to quantify whether its accuracy as a stock discriminant between genetically distinct populations can fluctuate as a response to extrinsic drivers. Further work should now investigate if the accuracy of the discrimination of known stocks also remains stable over time. Such study would provide interesting conclusions on the reliability over time of stock discrimination using otolith shape, and consequently on how different local stocks might respond to changes in the environment or in their population structure. However, even if these results are promising, some limitations must be considered. Because this study focused on within-stock stability of otolith shape, it comprised fish from a single stock unit and only individuals whose origin was previously determined with certainty based on inner otolith morphological features. While the separation criteria used by age readers have shown high agreement with independent genetic analyses (Berg et al., 2005), it is possible that excluding the individuals of uncertain origin in our analysis consequently removed some of the variability occurring within the NEA cod stock. However, mean annual percentage of uncertain NEA cod in the archives was low at only 5.24% (see Supplementary material 1). During the 1980s, the proportion of cod labelled as uncertain NEA was higher, primarily due to the changes in the readers’ organization. Besides, the discrimination criterion discussed here relies on the use of inner ring morphology, which is not directly comparable with the outer contour used in our analysis. Therefore, we presume that the potential effect of our sample selection is rather negligible when quantifying within-stock NEA cod otolith shape variability. Moreover, otoliths are physical structures whose features are defined in three dimensions, but the elliptical Fourier decomposition used here considers only the 2-dimensional projection of the otolith contour and in turn ignores potential differences in the three-dimensional shape. This method has been shown to downweigh the importance of higher-order harmonics, having most of the shape variance expressed in the first harmonics (Harbitz and Albert, 2015). This essentially means that many studies are biased against higher-order variations, which could mask subtle shape differences and factors of biological importance. In addition, elliptical Fourier descriptors are not completely independent from one another since each subsequent harmonic is built on the previous one, meaning that inherent biases could confound the statistical evaluation of shape differences (Haines and Crampton, 2000). While it remains one of the most widespread and accessible tools for shape analysis, EFA is ultimately one technique among many. For example, both geometric morphometrics and wavelet analysis show promising results and solve some of the Fourier decomposition shortcomings (Monteiro et al., 2005; Tuset et al., 2006; Vignon and Morat, 2010; Sadighzadeh et al., 2014). Results presented herein are thus constrained to the EFA approach applied on the selected stock and might not be applicable on a more general scale without further research. Sources of shape variation In accordance with previous studies, individual-specific variables contributed most to the observed variability in individual otolith shape (Campana and Casselman, 1993; Hüssy, 2008; Capoccioni et al., 2011). Body length, in particular, had a significant effect on shape and improved the models substantially, which is why most of the variation of otolith shape is usually attributed to allometry (Simoneau et al., 2000; Monteiro et al., 2005). By using fish of a single age-class throughout the analysis, differences in body length served as a proxy for different growth rates. These findings thus confirm that, while the ontogenetic trajectory of otolith shape is consistent between conspecific adults, changes in rate or timing of fish growth can lead to significant inter-individual differences in otolith morphology. Otolith orientation had also a significant effect on otolith shape. If bilateral asymmetry is prevalent, results of shape analyses can be affected by using different combinations of left, right, or both orientations (Mahé et al., 2019). Our findings could thus support the use of only left or right otoliths in further analyses. However, due to the age estimation routines, only one whole otolith from each fish was available and used in our analysis. The proportion of variability explained by otolith orientation was consequently inflated, since it also reflected inter-individual differences in shape that were not explained by the other predictors included in our models. In contrast, most extrinsic variables had a significant albeit much weaker effect on shape, which could indicate that changes in the environment or population dynamics over time account for less than the inter-individual variability occurring within a population. Although the reconstructed average cohort shapes were closely similar, their analysis revealed more clearly the effect of the different environmental variables. Interestingly, it also highlighted a large joint contribution of environmental variables and average cohort fish length, which suggests that environmental factors could exert an indirect influence on otolith shape through changes in fish growth rates. The influence of temperature on fish growth has been well documented (Brander, 2000; Björnsson et al., 2001) and has been linked to variations in otolith shape (Lombarte and Lleonart, 1993; Cardinale et al., 2004). However, because somatic growth is also consumption-dependent, feeding has been proposed as another important source of variability in otolith shape (for cod in Cardinale et al., 2004; Hüssy, 2008; for various flatfishes and roundfishes in Mille et al., 2016). Food availability depends in part on environmental conditions regulating the abundance of prey populations, as well as on competition for the access to food resources. Consequently, the relatively higher effect of total and spawning cod biomass on otolith shape, as well as the large joint contribution of extrinsic and intrinsic factors observed on average cohort shapes could, in fact, be representative of changes in growth rate driven by density-dependent effects on prey availability. However, the influence of such effects has long been debated for Gadus morhua due to the numerous other factors influencing food availability, which were not tested for here (Jørgensen, 1992; Brander, 2007). Further research should focus on the extent to which density-dependence can affect cod growth and in turn otolith growth. Ultimately, these extrinsic factors only explained a minor portion of the observed variability and the otolith shape of NEA cod appeared mostly stable within the stock. However, it should be stressed that this study focused on long-term, between-cohort changes in otolith shape in relation to the environment. A recent study has challenged the common paradigm that only long-term environmental differences drive changes in otolith shape by showing how acute, short-term episodes of environmental disturbances during early life can have similar results (Vignon, 2018). As such, significant otolith shape differences could likely occur within cohorts of the same stock following extreme climatic episodes. Further work is thus needed to fully understand the mechanisms behind changes in otolith shape in contrasted environmental conditions and at different scales, both temporally and spatially. Exploration of optimal environmental signals Identifying the critical window of action for an environmental factor is essential to interpreting biological responses of interest, where simplistic or wrongly identified signals can often depart from biological realities (van de Pol et al., 2016). In that context, the sliding window analysis provides a promising tool to investigate and identify optimal signals with a systematic and statistically sound approach. However, it is essential to interpret the results with caution considering the exploratory nature of technique. Constrained ordination methods like redundancy analysis do not provide a likelihood statistic and it is not possible to calculate an associated Akaike information criterion (AIC) or deviance. As explained by Oksanen et al. (2019), the “deviance” parameter generated for multivariate analysis used in this study is, in fact, akin to the residual sum of squares. According to both Borcard et al. (2018) and Oksanen et al. (2019), this criterion is not completely trustworthy and its interpretation can be quite liberal. The authors consequently stressed that it should be used carefully, and other criteria should be considered when available. It is nonetheless worth noting that this “deviance” parameter is also often used for the stepwise reduction of variables commonly seen in multivariate models, where similar caveats apply. In the absence of a likelihood statistic suitable for model comparison within this sliding window approach (Bailey and van de Pol, 2016; van de Pol et al., 2016), this artificial deviance criterion still provided a convenient tool for identifying optimal models. Besides, the sliding window analysis carried out in this study investigated factors that ultimately explained very low variability in the otolith shape decomposed into Fourier descriptors. Therefore, the models generated during the sliding window analysis identified optimum signals based on a limited proportion of the overall variance, which is statistically challenging. In this context, biological interpretations can be difficult and conclusions should be drawn carefully. Prior to analysis, it was hypothesized that environmental variables associated with the stock dynamics would have the most prominent effects during the early life of fish, since it is then most sensitive to cannibalism and lasting effects of early competition on its ontogeny (Folkvord and Otterå, 1993; Yaragina et al., 2009). However, total stock and spawning biomass critical windows were both identified for a period beginning at the transition from juvenile stage to adult and ending around the age cod reaches sexual maturity (Jørgensen, 1990). Because total stock and spawning biomass estimations do not comprise the youngest fish that have not been recruited yet, both variables could be more significant at later stages because they constitute proxies of the magnitude of the competition observed further during the fish life cycle. The density-dependent effects in the early life of fish are, therefore, not captured properly and identification of these relationships may be impossible. However, multiple studies have shown a strong ontogenetic development of Atlantic cod diet, in terms of both size and age. Link and Garrison (2002) observed ontogenetic shifts in Northeast Atlantic cod diet, wherein juveniles fed on small pelagic invertebrates, medium cod on a mix of benthic invertebrates and fish, and bigger cod on a larger amount of fish. Similarly, Jaworski and Ragnarsson (2006) described comparable shifts in the diet of Atlantic cod around Iceland. The present results could indicate that the density-dependent effects on fish growth (and consequently on otolith growth and shape) may be more significant during early adult life up to sexual maturity, when growth will then slow down. Moreover, randomization tests for both predictors confirmed that the identified optimal signal was strongly supported, which further justifies its biological relevance. In contrast, the optimal window identified for temperature revealed a narrower, earlier period of influence beginning prior to hatching and covering the juvenile stage of an individual’s life. Temperature is known as one of the major drivers of biological changes in fish such as growth rate and maturity (Brander, 2000), and its influence on otolith growth can be both indirect through changes of fish growth or through its direct effects on material deposition (Campana and Casselman, 1993; Campana, 1999). Indeed, otolith growth is an acellular process under different regulations than somatic growth (Simkiss, 1974), where biomineralization depends in part on the local chemistry (Borelli et al., 2003). Because these processes are directly affected by temperature (Casselman, 1990; Lombarte and Lleonart, 1993) changes in environmental conditions during early life, when both fish and otolith growth are maximal, could have a more significant influence on the final morphology of adult fish otoliths. Alternatively, early life temperature exposure could also have an indirect influence on fish and otolith growth through its effects on prey availability. Because Atlantic cod diet at the youngest life stages mostly consists of zooplankton for the larvae and small invertebrates for the juveniles (Link and Garrison, 2002; Jaworski and Ragnarsson, 2006), changes in temperature during these periods could lead to a mismatch with prey abundance and result in higher competition and poorer growth (Cushing, 1990; Rogers et al., 2011). Indeed, Hüssy (2008) showed that feeding level directly influenced otolith shape of young-of-year cod by affecting the number and dimensions of growth centres and lobes forming the crenulated edges of cod otolith. Differences in growth rate and otolith formation linked to temperature and feeding during early life could, therefore, have long-lasting effects on the shape of otoliths (Irgens et al., 2017) between individuals of the same stock. Interestingly, similar results were experimentally found by Cardinale et al. (2004). While this temperature signal was significant it should, however, be more carefully interpreted, because randomization tests showed a higher probability of getting a model supported by chance than for the previously identified optimal windows. Likewise, the average cohort shape analysis showed a high probability of the temperature signal to be a product of chance. However, this window was closely similar to the one previously identified in the analysis of individual shapes, and the signal for both total stock and SSB remained strongly supported by randomization tests. Furthermore, the reconstructed average cohort shapes were shown to be almost identical, so it can be hypothesized that the differences in explained variance between different tested time windows were, in fact, simply too low to be accurately separated. Because they closely resembled the windows identified in the analysis of individual shapes, including the temperature window in the models can still be justified since it significantly improved them. The identification of optimal windows using average cohort shapes should, however, be carefully considered, as it may have a limited statistical power in comparison with the individual analysis due to their low inter-variability. Despite its exploratory nature, the sliding window analysis offers compelling evidence that models accounting for environmental factors should be contrasted to best represent biological realities. The present study shows the potential of this approach not only to investigate environmental predictors selected a priori based on formulated hypotheses, but also to uncover potentially unknown signals of biological relevance. These results encourage future research to delve deeper into newer techniques that challenge common assumptions, especially when investigating the biological responses to environmental changes. Conclusions The present study generated a valuable continuous database of otolith shape extending over 82 years (1933–2015), the longest to date for this specific stock of Atlantic cod. This first long-term analysis of NEA cod otolith shape stability revealed that the general within-stock morphology remains largely unchanged through time and changing conditions. Similar work should now be conducted to identify shape stability within other species or stocks where shape analysis is a critical asset to fisheries management. In addition, the hierarchical modelling approach helped disentangle and partition the effects of different suspected sources of variation, further quantifying the resilience of otolith shape to changing environmental conditions. The sliding window approach and the effects of identified optimal signals on shape give new insights into the mechanisms behind otolith shape variability, which further contributes towards a consensus on the biological interpretation of shape differences. Supplementary data Supplementary material is available at the ICESJMS online version of the manuscript. Acknowledgements We want to thank Erlend Langhelle and Åse Husebø from the Institute of Marine Research (Norway) for their assistance with locating and imaging the otoliths. We benefitted from valuable comments and discussions during the 2019 International Sclerochronology Conference in Croatia. The authors are also grateful to the anonymous reviewers for their comments and suggestions which helped to improve this manuscript. Funding Funding for this work was provided by the Icelandic Research Fund Grant 173906-051 to S.E. Campana. Author contributions C.D., A.J.G., and S.S. conceived and designed this study. C.D. collected and processed the samples. C.D. and S.S. conducted data analysis and interpretation. C.D. wrote the manuscript. A.J.G, S.S., and J.A.G. edited and revised the manuscript. References Bailey L. D. , van de Pol M. 2016 . climwin: an R toolbox for climate window analysis . PLoS One , 11 : e0167980 . Google Scholar Crossref Search ADS PubMed WorldCat Barthelme S. 2019 . imager: Image Processing Library Based on ‘CImg’. R Package Version 0.41.2. https://CRAN.R-project.org/package=imager (last accessed 1 September 2019). Begg G. A. , Brown R. W. 2000 . Stock identification of haddock Melanogrammus aeglefinus on Georges Bank based on otolith shape analysis . Transactions of the American Fisheries Society , 129 : 935 – 945 . Google Scholar Crossref Search ADS WorldCat Berg E. , Sarvas T. H. , Harbitz A. , Fevolden S. E. , Salberg A. B. 2005 . Accuracy and precision in stock separation of north-east Arctic and Norwegian coastal cod by otoliths—comparing readings, image analyses and a genetic method . Marine and Freshwater Research , 56 : 753 – 762 . Google Scholar Crossref Search ADS WorldCat Björnsson B. , Steinarsson A. , Oddgeirsson M. 2001 . Optimal temperature for growth and feed conversion of immature cod (Gadus morhua L.) . ICES Journal of Marine Science , 58 : 29 – 38 . Google Scholar Crossref Search ADS WorldCat Bochkov Y. A. 1982 . Water temperature in the 0-200 m layer in the Kola-Meridian in the Barents Sea, 1900-1981 . Sb. Nauchn. Trud. PINRO , 46 : 113 – 122 . OpenURL Placeholder Text WorldCat Bonhomme V. , Picq S. , Gaucherel C. , Claude J. 2014 . Momocs: outline analysis using R . Journal of Statistical Software , 56 : 1 – 24 . Google Scholar Crossref Search ADS WorldCat Borcard D. , Gillet F. , Legendre P. 2018 . Numerical Ecology with R . Springer Science, New York. 444 pp. Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Borelli G. , Mayer-Gostan N. , Merle P. L. , De Pontual H. , Boeuf G. , Allemand D. , Payan P. 2003 . Composition of biomineral organic matrices with special emphasis on turbot (Psetta maxima) otolith and endolymph . Calcified Tissue International , 72 : 717 – 725 . Google Scholar Crossref Search ADS PubMed WorldCat Brander K. 2000 . Effects of environmental variability on growth and recruitment in cod (Gadus morhua) using a comparative approach . Oceanologica Acta , 23 : 485 – 496 . Google Scholar Crossref Search ADS WorldCat Brander K. M. 2007 . The role of growth changes in the decline and recovery of North Atlantic cod stocks since 1970 . ICES Journal of Marine Science , 64 : 211 – 217 . Google Scholar Crossref Search ADS WorldCat Campana S. E. 1999 . Chemistry and composition of fish otoliths: pathways, mechanisms and applications . Marine Ecology Progress Series , 188 : 263 – 297 . Google Scholar Crossref Search ADS WorldCat Campana S. E. 2001 . Accuracy, precision and quality control in age determination, including a review of the use and abuse of age validation methods . Journal of Fish Biology , 59 : 197 – 242 . Google Scholar Crossref Search ADS WorldCat Campana S. E. 2005 . Otolith science entering the 21st century . Marine and Freshwater Research , 56 : 485 – 495 . Google Scholar Crossref Search ADS WorldCat Campana S. E. , Casselman J. M. 1993 . Stock discrimination using otolith shape analysis . Canadian Journal of Fisheries and Aquatic Sciences , 50 : 1062 – 1083 . Google Scholar Crossref Search ADS WorldCat Capoccioni F. , Costa C. , Aguzzi J. , Menesatti P. , Lombarte A. , Ciccotti E. 2011 . Ontogenetic and environmental effects on otolith shape variability in three Mediterranean European eel (Anguilla anguilla, L.) local stocks . Journal of Experimental Marine Biology and Ecology , 397 : 1 – 7 . Google Scholar Crossref Search ADS WorldCat Cardinale M. , Doering-Arjes P. , Kastowsky M. , Mosegaard H. 2004 . Effects of sex, stock, and environment on the shape of known-age Atlantic cod (Gadus morhua) otoliths . Canadian Journal of Fisheries and Aquatic Sciences , 61 : 158 – 167 . Google Scholar Crossref Search ADS WorldCat Casselman J. M. 1990 . Growth and relative size of calcified structures of fish . Transactions of the American Fisheries Society , 119 : 673 – 688 . Google Scholar Crossref Search ADS WorldCat Castonguay M. , Simard P. , Gagnon P. 1991 . Usefulness of Fourier analysis of otolith shape for Atlantic Mackerel (Scomber scombrus) stock discrimination . Canadian Journal of Fisheries and Aquatic Sciences , 48 : 296 – 302 . Google Scholar Crossref Search ADS WorldCat Cook B. I. , Wolkovich E. M. , Parmesan C. 2012 . Divergent responses to spring and winter warming drive community level flowering trends . Proceedings of the National Academy of Sciences , 109 : 9000 – 9005 . Google Scholar Crossref Search ADS WorldCat Cushing D. H. 1990 . Plankton production and year-class strength in fish populations: an update of the match/mismatch hypothesis . Advances in Marine Biology , 26 : 249 – 293 . Google Scholar Crossref Search ADS WorldCat Dahle G. , Johansen T. , Westgaard J.-I. , Aglen A. , Glover K. A. 2018 . Genetic management of mixed-stock fisheries “real-time”: the case of the largest remaining cod fishery operating in the Atlantic in 2007–2017 . Fisheries Research , 205 : 77 – 85 . Google Scholar Crossref Search ADS WorldCat Dippner J. W. , Ottersen G. 2001 . Cod and climate variability in the Barents Sea . Climate Research , 17 : 73 – 82 . Google Scholar Crossref Search ADS WorldCat Folkvord A. , Otterå H. 1993 . Effects of initial size distribution, day length, and feeding frequency on growth, survival, and cannibalism in juvenile Atlantic cod (Gadus morhua L.) . Aquaculture , 114 : 243 – 260 . Google Scholar Crossref Search ADS WorldCat Friedland K. D. , Reddin D. G. 1994 . Use of otolith morphology in stock discriminations of Atlantic salmon (Salmo salar) . Canadian Journal of Fisheries and Aquatic Sciences , 51 : 91 – 98 . Google Scholar Crossref Search ADS WorldCat Haines A. J. , Crampton J. S. 2000 . Improvements to the method of Fourier shape analysis as applied in morphometric studies . Palaeontology , 43 : 765 – 783 . Google Scholar Crossref Search ADS WorldCat Hammer C. , Zimmermann C. 2005 . Chapter 30—The role of stock identification in formulating fishery management advice . In Stock Identification Methods , pp. 631–658. Ed. by S. X. Cadrin, K. D. Friedland and J. R. Waldman. Elsevier Academic Press, Burlington, USA. OpenURL Placeholder Text WorldCat Harbitz A. , Albert O. T. 2015 . Pitfalls in stock discrimination by shape analysis of otolith contours . ICES Journal of Marine Science , 72 : 2090 – 2097 . Google Scholar Crossref Search ADS WorldCat Hüssy K. 2008 . Otolith shape in juvenile cod (Gadus morhua): ontogenetic and environmental effects . Journal of Experimental Marine Biology and Ecology , 364 : 35 – 41 . Google Scholar Crossref Search ADS WorldCat Hylen A. 2002 . Fluctuations in abundance of Northeast Arctic cod during the 20th century . ICES Marine Science Symposia , 215 : 543 – 550 . OpenURL Placeholder Text WorldCat ICES. 2018 . Report of the Arctic Fisheries Working Group (AFWG), 18–24 April 2018, Ispra, Italy. ICES Document CM 2018/ACOM:06. Irgens C. , Kjesbu O. S. , Folkvord A. 2017 . Ontogenetic development of otolith shape during settlement of juvenile Barents Sea cod (Gadus morhua) . ICES Journal of Marine Science , 74 : 2389 – 2397 . Google Scholar Crossref Search ADS WorldCat Jaworski A. , Ragnarsson S. Á. 2006 . Feeding habits of demersal fish in Icelandic waters: a multivariate approach . ICES Journal of Marine Science , 63 : 1682 – 1694 . Google Scholar Crossref Search ADS WorldCat Johansen T. , Westgaard J.-I. , Seliussen B. B. , Nedreaas K. , Dahle G. , Glover K. A. , Kvalsund R. et al. 2018 . “Real-time” genetic monitoring of a commercial fishery on the doorstep of an MPA reveals unique insights into the interaction between coastal and migratory forms of the Atlantic cod . ICES Journal of Marine Science , 75 : 1093 – 1104 . Google Scholar Crossref Search ADS WorldCat Jónsdóttir I. G. , Campana S. E. , Marteinsdottir G. 2006 . Otolith shape and temporal stability of spawning groups of Icelandic cod (Gadus morhua L.) . ICES Journal of Marine Science , 63 : 1501 – 1512 . Google Scholar Crossref Search ADS WorldCat Jørgensen T. 1990 . Long-term changes in age at sexual maturity of Northeast Arctic cod (Gadus morhua L.) . ICES Journal of Marine Science , 46 : 235 – 248 . Google Scholar Crossref Search ADS WorldCat Jørgensen T. 1992 . Long-term changes in growth of Northeast Arctic cod (Gadus morhua) and some environmental influences . ICES Journal of Marine Science , 49 : 263 – 278 . Google Scholar Crossref Search ADS WorldCat Kruuk L. E. B. , Osmond H. L. , Cockburn A. 2015 . Contrasting effects of climate on juvenile body size in a Southern Hemisphere passerine bird . Global Change Biology , 21 : 2929 – 2941 . Google Scholar Crossref Search ADS PubMed WorldCat Kuhl F. P. , Giardina C. R. 1982 . Elliptic Fourier features of a closed contour . Computer Graphics and Image Processing , 18 : 236 – 258 . Google Scholar Crossref Search ADS WorldCat Legendre P. , Legendre L. F. 2012 . Numerical Ecology . Elsevier , Oxford. 281 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Lestrel P. E. 1997 . Fourier Descriptors and Their Applications in Biology . Cambridge University Press , Cambridge . 460 pp. Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Link J. S. , Garrison L. P. 2002 . Trophic ecology of Atlantic cod (Gadus morhua) on the northeast US continental shelf . Marine Ecology Progress Series , 227 : 109 – 123 . Google Scholar Crossref Search ADS WorldCat Lombarte A. , Lleonart J. 1993 . Otolith size changes related with body growth, habitat depth and temperature . Environmental Biology of Fishes , 37 : 297 – 306 . Google Scholar Crossref Search ADS WorldCat Mahé K. , Ider D. , Massaro A. , Hamed O. , Jurado-Ruzafa A. , Gonçalves P. , Anastasopoulou A. et al. 2019 . Directional bilateral asymmetry in otolith morphology may affect fish stock discrimination based on otolith shape analysis . ICES Journal of Marine Science , 76 : 232 – 243 . Google Scholar Crossref Search ADS WorldCat Mahe K. , Oudard C. , Mille T. , Keating J. , Gonçalves P. , Clausen L. W. , Petursdottir G. et al. 2016 . Identifying blue whiting (Micromesistius poutassou) stock structure in the Northeast Atlantic by otolith shape analysis . Canadian Journal of Fisheries and Aquatic Sciences , 73 : 1363 – 1371 . Google Scholar Crossref Search ADS WorldCat Mille T. , Mahé K. , Cachera M. , Villanueva M. C. , Pontual H. D. , Ernande B. 2016 . Diet is correlated with otolith shape in marine fish . Marine Ecology Progress Series , 555 : 167 – 184 . Google Scholar Crossref Search ADS WorldCat Monteiro L. R. , Beneditto A. P. M. D. , Guillermo L. H. , Rivera L. A. 2005 . Allometric changes and shape differentiation of sagitta otoliths in sciaenid fishes . Fisheries Research , 74 : 288 – 299 . Google Scholar Crossref Search ADS WorldCat Oksanen J. , Blanchet G. F. , Friendly M. , Kindt R. , Legendre P. , McGlinn D. , Minchin P. R. et al. 2019 . vegan: Community Ecology Package. R Package Version 2.5-4. https://CRAN.R-project.org/package=vegan (last accessed 1 September 2019). Peterson B. G. , Carl P. , Boudt K. , Bennett R. , Ulrich J. , Zivot E. et al. 2019 . Package ‘PerformanceAnalytics’. R Team Cooperation. http://www.r-forge.r-project.org/projects/returnanalytics (last accessed 1 September 2019). Petursdottir G. , Begg G. A. , Marteinsdottir G. 2006 . Discrimination between Icelandic cod (Gadus morhua L.) populations from adjacent spawning areas based on otolith growth and shape . Fisheries Research , 80 : 182 – 189 . Google Scholar Crossref Search ADS WorldCat R Core Team. 2016 . A Language and Environment for Statistical Computing . R Foundation for Statistical Computing , Vienna, Austria . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Roberts A. M. I. , Tansey C. , Smithers R. J. , Phillimore A. B. 2015 . Predicting a change in the order of spring phenology in temperate forests . Global Change Biology , 21 : 2603 – 2611 . Google Scholar Crossref Search ADS PubMed WorldCat Robichaud D. , Rose G. A. 2004 . Migratory behaviour and range in Atlantic cod: inference from a century of tagging . Fish and Fisheries , 5 : 185 – 214 . Google Scholar Crossref Search ADS WorldCat Rogers L. A. , Stige L. C. , Olsen E. M. , Knutsen H. , Chan K.-S. , Stenseth N. C. 2011 . Climate and population density drive changes in cod body size throughout a century on the Norwegian coast . Proceedings of the National Academy of Sciences , 108 : 1961 – 1966 . Google Scholar Crossref Search ADS WorldCat Rollefsen G. 1934 . The cod otolith as a guide to race, sexual development and mortality . Rapports et Proces-Verbaux des Réunions Conseil International pour l’Exploration de la Mer , 88 : 1 – 6 . OpenURL Placeholder Text WorldCat Sadighzadeh Z. , Valinassab T. , Vosugi G. , Motallebi A. A. , Fatemi M. R. , Lombarte A. , Tuset V. M. 2014 . Use of otolith shape for stock identification of John’s snapper, Lutjanus johnii (Pisces: Lutjanidae), from the Persian Gulf and the Oman Sea . Fisheries Research , 155 : 59 – 63 . Google Scholar Crossref Search ADS WorldCat Salvanes A. G. V. , Skjæraasen J. E. , Nilsen T. 2004 . Sub-populations of coastal cod with different behaviour and life-history strategies . Marine Ecology Progress Series , 267 : 241 – 251 . Google Scholar Crossref Search ADS WorldCat Simkiss K. , 1974. Calcium metabolism of fish in relation to ageing. In The Ageing of Fish, pp. 1–12. Ed. by T. B. Bagenal. Unwin Brothers, Old Woking. Simoneau M. , Casselman J. M. , Fortin R. 2000 . Determining the effect of negative allometry (length/height relationship) on variation in otolith shape in lake trout (Salvelinus namaycush), using Fourier-series analysis . Canadian Journal of Zoology , 78 : 1597 – 1603 . Google Scholar Crossref Search ADS WorldCat Smoliński S. 2019 . Sclerochronological approach for the identification of herring growth drivers in the Baltic Sea . Ecological Indicators , 101 : 420 – 431 . Google Scholar Crossref Search ADS WorldCat Stransky C. , Baumann H. , Fevolden S.-E. , Harbitz A. , Høie H. , Nedreaas K. H. , Salberg A.-B. et al. 2008 . Separation of Norwegian coastal cod and Northeast Arctic cod by outer otolith shape analysis . Fisheries Research , 90 : 26 – 35 . Google Scholar Crossref Search ADS WorldCat Stransky C. , MacLellan S. E. 2005 . Species separation and zoogeography of redfish and rockfish (genus Sebastes) by otolith shape analysis . Canadian Journal of Fisheries and Aquatic Sciences , 62 : 2265 – 2276 . Google Scholar Crossref Search ADS WorldCat Tereshchenko V. V. 1996 . Seasonal and Year-to-year Variations of Temperature and Salinity along the Kola Meridian Transect. ICES Document CM 1996/C: 11. 24 pp. Tracey S. R. , Lyle J. M. , Duhamel G. 2006 . Application of elliptical Fourier analysis of otolith form as a tool for stock identification . Fisheries Research , 77 : 138 – 147 . Google Scholar Crossref Search ADS WorldCat Tuset V. M. , Rosin P. L. , Lombarte A. 2006 . Sagittal otolith shape used in the identification of fishes of the genus Serranus . Fisheries Research , 81 : 316 – 325 . Google Scholar Crossref Search ADS WorldCat van de Pol M. , Bailey L. D. , McLean N. , Rijsdijk L. , Lawson C. R. , Brouwer L. 2016 . Identifying the best climatic predictors in ecology and evolution . Methods in Ecology and Evolution , 7 : 1246 – 1257 . Google Scholar Crossref Search ADS WorldCat Vignon M. 2015 . Disentangling and quantifying sources of otolith shape variation across multiple scales using a new hierarchical partitioning approach . Marine Ecology Progress Series , 534 : 163 – 177 . Google Scholar Crossref Search ADS WorldCat Vignon M. 2017 . Investigating morphospace occupation in multi-scale ecological and evolutionary data using regression tree: case studies and perspectives . Evolutionary Biology , 44 : 120 – 134 . Google Scholar Crossref Search ADS WorldCat Vignon M. 2018 . Short-term stress for long-lasting otolith morphology—brief embryological stress disturbance can reorient otolith ontogenetic trajectory . Canadian Journal of Fisheries and Aquatic Sciences , 75 : 1713 – 1722 . Google Scholar Crossref Search ADS WorldCat Vignon M. , Morat F. 2010 . Environmental and genetic determinant of otolith shape revealed by a non-indigenous tropical fish . Marine Ecology Progress Series , 411 : 231 – 241 . Google Scholar Crossref Search ADS WorldCat Wickham H. 2017 . tidyverse: Easily Install and Load the ‘tidyverse’. R Package Version 1.2.1. https://CRAN.R-project.org/package=tidyverse (last accessed 1 September 2019). Yaragina N. A. , Bogstad B. , Kovalev Y. A. 2009 . Variability in cannibalism in Northeast Arctic cod (Gadus morhua) during the period 1947–2006 . Marine Biology Research , 5 : 75 – 85 . Google Scholar Crossref Search ADS WorldCat © International Council for the Exploration of the Sea 2020. All rights reserved. For permissions, please email: [email protected] This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Parsa,, Mahdi;Emery, Timothy, J;Williams, Ashley, J;Nicol,, Simon

doi: 10.1093/icesjms/fsaa020pmid: N/A

Abstract Minimizing fishing-induced mortality on bycatch and endangered, threatened or protected species is a necessity for fisheries managers. Estimating individual vessel bycatch rates by dividing the amount of bycatch by effort (nominal rate) can be biased, as it does not consider effort heterogeneity within the fleet and ignores prior knowledge of fleet bycatch rates. We develop an empirical Bayesian approach for estimating individual vessel and fleet bycatch rates that: (i) considers effort heterogeneity among vessels and; (ii) pools data from similar vessels for more accurate estimation. The proposed standardized bycatch rate of a vessel is, therefore, the weighted average of the pool rate and nominal rate of the vessel; where the weights are functions of the vessel’s fishing effort and a constant estimated from the model. We apply this inference method to the estimation of seabird bycatch rates in the component of the Australian Eastern Tuna and Billfish Fishery targeting yellowfin tuna. We illustrate the capability of the method for providing fishery managers with insights on fleet-wide bycatch mitigation performance and the identification of outperforming and underperforming vessels. This method can also be used by fishery managers to develop fleet-wide performance measures or quantitative evaluation standards. Introduction Global fisheries bycatch in wild-capture fisheries is an issue of growing concern (Diamond, 2004; Gilman et al., 2008). Species that have little or no economic value to fishers (e.g. due to their small size); prohibited species (e.g. those managed in other fisheries); regulatory discards (e.g. species below or above the size limit); or endangered, threatened or protected (ETP) species (e.g. marine turtles, seabirds) are all examples of bycatch species (Diamond, 2004). For this article, we refer hereafter to bycatch species as those species that are caught and subsequently discarded at sea, or in the case of ETP species, interacted with at sea. While the 1982 United Nations Convention of the Law of the Sea under Article 61 requires signatories to determine the biological and ecological impacts of fishing on non-target (bycatch) species, this can be difficult for most commercial fisheries that lack fishery-dependent data. As reported by Tuck (2011), bycatch data are often limited due to inadequate and incomplete information on vessel characteristics, fishing effort, and species composition. Many species are under- or over-reported, non-reported, or misreported in fishery logbooks (Walsh et al., 2002; Walsh et al., 2005; Sampson, 2011; Mangi et al., 2016; Macbeth et al., 2018). For example, in an examination of catch rates for blue shark (Prionace glauca), Walsh et al. (2002) found that underreported catches in fishery logbooks were due to fishers being too busy to report incidental catches. In a similar study examining the catch rates for blue marlin (Makaira nigricans), Walsh et al. (2005) observed that fishers tended to over-report catches due to misidentifying striped marlin (Tetrapturus audax) and shortbill spearfish (Tetrapturus angustirostris) as blue marlin. The inadequacies of fishery logbook data have often led decision-makers to use at-sea observer data as an alternative to quantify bycatch taken by commercial fisheries. However, at-sea observer data have its own suite of biases (Benoît and Allard, 2009; Faunce and Barbeaux, 2011; Wakefield et al., 2018) and any extrapolations of at-sea observer data at low levels of coverage are likely to produce imprecise and inaccurate results when capture of a species is a rare occurrence (Wakefield et al., 2018). Despite the issues associated with logbook data, it often remains the principal source of information on fishery catch and effort due to many management authorities requiring vessels to fill out their logbook as a condition of their licence or permit (Sampson, 2011). Access to fishery logbook data allows the nominal discard rate for bycatch species to be calculated at an individual vessel or fleet level. This is often done by dividing the amount of bycatch by the total effort for a given vessel. This is termed the “nominal” estimate. This vessel-level estimation could be unbiased if there are sufficient observations (i.e. adequate sample size), and fishers have not changed their fishing practices over the time period assessed. However, this is often not the case, as different vessels enter and exit the fishery through time and change their fishing practices, influencing catchability (Tuck, 2011). Furthermore, consider two longline vessels with the same standard seabird bycatch rate of zero (0.0 bycatch per 1000 hooks), where vessel 1 expended a significantly greater amount of effort compared with vessel 2. Calculation of the nominal estimate would suggest that both vessels are performing identically; however, from the perspective of a fishery manager, vessel 1 is outperforming vessel 2 since there has been no bycatch recorded with a substantially greater exposure to risk (i.e. effort). Moreover, a fishery manager is more confident in the bycatch rate of vessel 1, simply due to the greater level of effort expended compared with vessel 2, whose zero-bycatch rate could simply be due to chance through limited exposure. The nominal estimate also only uses each vessel’s information for estimating the rate and ignores other available information (e.g. effort data) from “similar” vessels in each fleet or fishery. Given these limitations, we propose a “standardized” estimate using an empirical Bayesian approach that considers effort heterogeneity among the fleet and pools data from “similar” vessels for rate estimation. Similar vessels are defined as those that share comparable fishing behaviour patterns [e.g. “fishing styles” after Boonstra and Hentati-Sundberg (2016) or “fishing tactics” after Pelletier and Ferraris (2000)] and can be pre-determined using variable quantitative or semi-quantitative methods based on the data from the commercial fishery or expert judgement, respectively. Vessel-, fleet- and fishery-level estimations of bycatch rates are sources of information that assist fisheries managers with monitoring the performance of bycatch mitigation measures. Vessel-level estimation may provide insight (through a targeted investigation) on why a vessel is underperforming (higher bycatch rate) or outperforming (lower bycatch rate) the fleet average (e.g. due to fishing in an area with the high abundance of protected species or appropriately deploying mitigation devices, respectively). Comparing the vessel-level estimated bycatch rates to the fleet-level estimate ensures that individual vessels are accountable for their actions and allows managers to set quantifiable bycatch thresholds for the fishery. Quantifiable measures, standards or reference points that guide expected levels of performance can create incentives for industry to reduce their bycatch rates through, for example altering fishing behaviour or adopting alternative bycatch mitigation technology (Diamond, 2004; Grafton et al., 2007; Kirby and Ward, 2014; Lent and Squires, 2017). When these performance standards create market-based incentives or disincentives (carrots and sticks) for industry, they have the potential to further improve fleet bycatch performance and reduce regulatory costs (Gjertsen et al., 2010; Pascoe et al., 2010). For example, in Australia, there is a Threat Abatement Plan (TAP) for seabirds, which sets a maximum permissible bycatch rate of 0.01 or 0.05 birds per 1000 hooks in various Australian Commonwealth fisheries (Commonwealth of Australia, 2018). Attached to this performance measure are criteria developed to guide the management response when the bycatch rate is exceeded, which may target individual vessels or the fleet and may have immediate economic costs (Commonwealth of Australia, 2018). In this article, we outline an inference method for calculating a model-estimated (standardized) bycatch rate for each vessel, which is the weighted average of the pool (fleet) rate and the nominal estimation rate of the individual vessel. Using an empirical Bayesian approach for the analysis of rare-event data is not new (Myers et al., 2002; Quigley et al., 2011) and has been shown to produce less biased and more consistent estimates of the probabilities of rare events compared with conventional statistical methods (Khakzad et al., 2014). We apply this method to a case study of seabird bycatch rates in the yellowfin tuna component of the Australian Eastern Tuna and Billfish Fishery (ETBF). We use the Australian ETBF as an example because we are confident that the fishery logbook data are the accurate representation of catch composition and bycatch of protected species in the years subsequent to the introduction of electronic monitoring technologies (Emery et al., 2019a). The results of the analysis are discussed in the context of (i) developing quantitative performance standards for bycatch species; (ii) reducing the transaction costs of management decision-making through a risk-based approach; and (iii) making fishers individually accountabile for their bycatch rates. Methodology Poisson–gamma model to estimate bycatch rates In our model, we assume that the amount of bycatch is approximately proportional to the total units of effort. This assumption is valid and is supported by the existing literature (Hatch, 2018) and the results of our study (see below). To estimate the standardized (seabird bycatch) rate of individual vessels, we develop a Poisson–gamma (Carlin and Louis, 2009) model considering two sources of uncertainties: (i) the uncertainties that arise from the lack of knowledge (e.g. the actual bycatch rate is not known), termed epistemic uncertainty, and (ii) uncertainty associated with natural variations in the sample (e.g. same amount of effort leads to a different amount of bycatch), termed aleatory uncertainties. Consequently, we use a gamma prior distribution to capture epistemic uncertainties within the pool of data to allow us to model the variation in true bycatch (actual seabird bycatch) rates, which are currently unknown. That is, we assume that the true bycatch rate of vessel i is a random variable with the gamma distribution of shape parameter α and scale parameter β ⁠. We denote it by λi∼gamma(α,β) ⁠, and the gamma probability density function can be expressed as the following equation. The mean of a gamma distribution is αβ ⁠, and here, we refer it as the pool rate. π(λi)=βαλiα-1e-βλiΓα,α>0,β>0,λi>0. (1) We later update the prior for each vessel to estimate the standardized bycatch rate. The updating process can be done quickly as the posterior of the gamma distribution remains in the gamma family, and we only need to update the shape and scale parameters. If we assume that n0 bycatch species were observed for E0 units of effort, Bayes’ theorem implies that the posterior distribution is of the form of the following equation: πλn0,E0=β+E0αλα+n0-1e-β+E0λΓα+n0,α, β,λ, E0>0,n0=0,1,2,3,…. (2) Assuming that the true bycatch rate Λi=λi for vessel i is constant for given Ei units of effort, we can then model the aleatory uncertainty in the bycatch rate through a Poisson probability distribution expressed in the following equation: PNi=niΛi=λi=λiEnie-λiEin!,Ei>0,λi>0,ni=0,1,2,…. (3) Since we do not know the true bycatch rate Λi for vessel i ⁠, we average the Poisson distributions, weighted against the prior distribution in the following equation: PNi=ni=∫0∞λiEinie-λiEini!βαλiα-1e-βλiΓαdλ,α>0, β>0, ni=0,1,2,…. (4) Greenwood and Yule (1920) proved that the distribution of Ni is Negative Binomial as shown in the following equation: PNi=ni=Γni+αΓαni!ββ+EiαEiβ+Eini,α>0, β>0, ni=0,1,2,…. (5) To estimate the parameters of the prior distribution, α,β ⁠, we use a genetic algorithm optimization method (implemented in MATLAB Global Optimization Toolbox) to maximize the natural logarithm of the marginal likelihood (LML) functions assuming that (pooled) data are generated from the Negative Binomial distribution of (5). Our choice of algorithm was informed by as follows: (i) there being no closed-form solution for finding maximum values of LML functions and (ii) the LML functions being highly nonlinear and nonconvex. Several methods have been proposed to construct a joint confidence region to address the uncertainty associated with the estimated prior parameters, such as the bootstrap method (Carlin and Gelfand, 1991), and using likelihood theory by assuming the negative of two times the natural logarithm of the relative marginal likelihood function has a chi-square distribution with two degrees of freedom (Basu and Rigdon, 1986). In this study, we used the second approach to construct a joint confidence interval for the maximum likelihood estimates and consequently the posterior mean (standardized) bycatch rate of each vessel. We let α^ and β^ are the estimated values of prior parameters and let vessel i interacts with ni bycatch species when Ei units of effort have been deployed. We estimate the standardized bycatch rate of vessel i ⁠, which is the posterior mean of λi as follows: Eλi|Ni=ni=∫0∞λiπ(λi|Ni=ni,α^,β^)dλi=α^+niβ^+Ei=α^β^(1-z)+niEiz, (6) where z=Eiβ^+Ei. The standardized bycatch rate can be interpreted as a weighted average of the pool (i.e. fleet) mean bycatch rate (⁠ α^/β^) and the nominal bycatch rate of the vessel (⁠ ni/Ei ⁠) where the weight is the function of a vessel’s fishing effort and a scale parameter of the posterior gamma distribution. Equation (6) also implies that when we have more experience (i.e. fishing effort) with a vessel (higher E ⁠), more weight will be allocated to the nominal rate, while for a vessel with less experience, more weight will be allocated to the pool rate. Application of the Poisson–gamma model to the Australian yellowfin tuna sub-fishery We apply this method to vessels in the yellowfin tuna sub-fishery of the Australian ETBF to illustrate how the method can provide fishery managers with insights on fleet-wide bycatch mitigation performance and identify non-performing vessels for targeted intervention. The ETBF is a pelagic longline fishery that operates within the Australian Exclusive Economic Zone and adjacent high sea waters targeting yellowfin tuna (Thunnus albacares), bigeye tuna (Thunnus obesus), albacore tuna (Thunnus alulunga), broadbill swordfish (Xiphias gladius), and striped marlin (T. audax). The ETBF operates from Cape York, east and south to the Victorian–South Australian border, including waters around Tasmania and the high seas of the Pacific Ocean (Figure 1a). In 2018, there were a total of 40 longline vessels active in the ETBF (Patterson et al., 2018). In the ETBF, vessels that have fished >30 days in the previous or current fishing season must have operational electronic monitoring technology installed. Figure 1. Open in new tabDownload slide Area and relative fishing intensity in the (a) eastern tuna and billfish fishery and (b) yellowfin tuna component of the eastern tuna and billfish fishery in 2016–2018 calendar years. Figure 1. Open in new tabDownload slide Area and relative fishing intensity in the (a) eastern tuna and billfish fishery and (b) yellowfin tuna component of the eastern tuna and billfish fishery in 2016–2018 calendar years. The yellowfin tuna sub-fishery of the Australian ETBF was differentiated from other sub-fisheries using a non-hierarchical clustering method, partitioning around medoids as similarly employed by Duarte et al. (2009) that identified structures within the data to quantitatively categorize individual fishing events to a particular métier (for more information on métier analysis, see Pelletier and Ferraris, 2000; Holley and Marchal, 2004). While the primary target species of the yellowfin tuna sub-fishery is yellowfin tuna, there is also a high proportion of oilfish (Ruvettus pretiosus) and striped marlin caught as by-products. The yellowfin tuna sub-fishery is a year-round fishery with most sets occurring between 7 and 9 a.m. off the New South Wales and Victorian State coastlines (Figure 1b). Typical gear characteristics include shallow setting with limited light stick use. In undertaking this analysis, we limit our study to the years 2016–2018 when electronic monitoring technologies were installed on all full-time ETBF vessels. This decision was based on recently published studies indicating that fishers have improved their logbook reporting of bycatch and protected species in these years, and there is high congruence between logbook and electronic monitoring analyst-reported seabird bycatch rates (Larcombe et al., 2016; Emery et al., 2019a, b). In 2016–2018, there were a total of 23, 29 and 26 longline vessels active, respectively, in this sub-fishery. Results Fishing effort in the yellowfin tuna sub-fishery There was high heterogeneity in the effort data for the 34 ETBF vessels operating in the yellowfin tuna sub-fishery during 2016–2018, with vessel_id 15 setting 216 000 hooks and vessel_id 6 and 21 just 1000 hooks, for example (Figure 2a). Furthermore, the amount of seabird bycatch varied among vessels with similar effort levels (Figure 2b). For example, vessel_id 16 and vessel_id 28 expended a similar amount of effort (160–180 000 hooks) in the yellowfin tuna sub-fishery between 2016 and 2018, but the number of recorded seabirds was different (six and one, respectively) (Figure 2b). Nevertheless, there was a positive linear correlation (Pearson’s r = 0.59, p = 0.00028) between the number of seabirds and the effort for each vessel. This result supports the assumption of proportionality between the amount of seabird bycatch and the amount of effort in the yellowfin tuna sub-fishery of the ETBF. Figure 2. Open in new tabDownload slide Total fishing effort (a) and amount of seabird bycatch (b) for a total of 34 vessels operating in the yellowfin tuna sub-fishery for the years 2016–2018. Figure 2. Open in new tabDownload slide Total fishing effort (a) and amount of seabird bycatch (b) for a total of 34 vessels operating in the yellowfin tuna sub-fishery for the years 2016–2018. Assessing seabird bycatch rates in the yellowfin tuna sub-fishery The mean seabird bycatch rate was 0.019 for the yellowfin tuna sub-fishery (i.e. average pool rate) based on (5), which was used in association with the nominal bycatch rate of the vessel in (6) to generate the standardized bycatch rate for each vessel. The standardized bycatch rate of a vessel with low levels of fishing effort was closer to the average pool rate, while the standardized bycatch rate of a vessel with high levels of fishing effort was closer to their nominal bycatch rate (Figure 3). Figure 3. Open in new tabDownload slide Standardized seabird bycatch rates for all 34 vessels in the yellowfin tuna sub-fishery for the years 2016–2018 plotted against their nominal bycatch rate. The size of each point represents the total effort of each vessel in ‘000s hooks. The red line is the identity line (1:1), and the blue line is the mean estimated bycatch rate for the fleet (i.e. average pool rate). Figure 3. Open in new tabDownload slide Standardized seabird bycatch rates for all 34 vessels in the yellowfin tuna sub-fishery for the years 2016–2018 plotted against their nominal bycatch rate. The size of each point represents the total effort of each vessel in ‘000s hooks. The red line is the identity line (1:1), and the blue line is the mean estimated bycatch rate for the fleet (i.e. average pool rate). The fit of the estimated predictive distribution model to the empirical data was robust (Figure 4). There was a good fit to the data in both the centre and right-hand tails of the distribution, while there was a slight overestimation and underestimation of the zero and one occurrences, respectively, on the left-hand tail of the distribution (Figure 4). The good fit to the upper right-hand tail of the distribution is very important since this has greater consequences for seabird populations if the true bycatch rate of a vessel is relatively high. Figure 4. Open in new tabDownload slide Hanging rootogram of the Poisson–gamma model fitted to seabird bycatch data for all 34 vessels in the yellowfin tuna sub-fishery for the years 2016–2018. The red line shows the expected amount of seabird bycatch estimated by the model, while the observed amount of seabird bycatch is shown as bars hanging from the red lines. The x-axis shows bins representing the nominal amount of seabird bycatch, while the y-axis shows the square root of the expected or observed amount of seabird bycatch. When the bar does not touch the x-axis (e.g. zero occurrences), it means that the amount of bycatch predicted by the model is higher than in the empirical data, while when the bar does touch the y-axis (e.g. one occurrence), it means that the amount of bycatch predicted by the model is lower than in the empirical data. Figure 4. Open in new tabDownload slide Hanging rootogram of the Poisson–gamma model fitted to seabird bycatch data for all 34 vessels in the yellowfin tuna sub-fishery for the years 2016–2018. The red line shows the expected amount of seabird bycatch estimated by the model, while the observed amount of seabird bycatch is shown as bars hanging from the red lines. The x-axis shows bins representing the nominal amount of seabird bycatch, while the y-axis shows the square root of the expected or observed amount of seabird bycatch. When the bar does not touch the x-axis (e.g. zero occurrences), it means that the amount of bycatch predicted by the model is higher than in the empirical data, while when the bar does touch the y-axis (e.g. one occurrence), it means that the amount of bycatch predicted by the model is lower than in the empirical data. It is evident that between 2016 and 2018 the average pool rate (red line in Figure 5) in the yellowfin tuna sub-fishery was below the maximum permissible bycatch rate of 0.05 seabird per 1000 hooks (blue line) recommended in the Australian Seabird TAP (Commonwealth of Australia, 2018) (Figure 5). However, there was a large variation among the 34 individual vessels, with some vessels having high standardized bycatch rates above the TAP (e.g. vessel_id 20, 22, and 32) and others having lower standardized bycatch rates (e.g. vessel_id 2, 5 and 8). The level of uncertainty in the estimated bycatch rates also varied substantially at the individual vessel level (Figure 5). Figure 5. Open in new tabDownload slide Standardized seabird bycatch rates for the 34 vessels in the yellowfin tuna sub-fishery for the years 2016–2018. The blue line represents the TAP recommended reference point (0.05 seabirds per 1000 hooks), and the red line represents the average pool rate. The grey shaded area represents the confidence interval for the estimated average pool rate. Figure 5. Open in new tabDownload slide Standardized seabird bycatch rates for the 34 vessels in the yellowfin tuna sub-fishery for the years 2016–2018. The blue line represents the TAP recommended reference point (0.05 seabirds per 1000 hooks), and the red line represents the average pool rate. The grey shaded area represents the confidence interval for the estimated average pool rate. Discussion Attaining robust estimates of bycatch rates in fisheries is a significant challenge due to their low (often rare in the case of ETP species) frequency of occurrence, leading to uncertainty in rate estimation, which can be a significant barrier to the development of effective mitigation strategies (Komoroske and Lewison, 2015; Martin et al., 2015; Suuronen and Gilman, 2019). Despite these challenges, fisheries managers are often required to make inferences about bycatch rates to inform their decision-making. This can lead to biased, imprecise estimates when using nominal estimation (dividing the total amount of bycatch by total effort) to determine the rate (Martin et al., 2015). By considering effort heterogeneity among vessels and pooling the data from homogenous vessels (vessels that share comparable fishing behavioural patterns), our model-estimated (standardized) bycatch rate overcomes some of the shortcomings of nominal estimation (Bishop et al., 2008). It also requires minimal data: only the total effort and amount of bycatch for each homogenous vessel within the timeframe of interest. This makes it more accessible to use in data-limited fisheries and easier for decision-makers to update and review regularly. Furthermore, by using Bayesian methods, which are well suited to the analysis of rare-event bycatch data, we can more fully integrate uncertainty, produce less volatile bycatch rate estimates, and enable evaluation of these estimates relative to existing performance measures (Gardner et al., 2008; Martin et al., 2015). We should emphasize that while other factors contribute to the bycatch rate, such as climate, location, food availability, and seasonality (Martin et al., 2015; Cortés et al., 2017), they were not considered in our model to ensure simplicity but could be incorporated as covariates in future modifications of this approach. Moreover, while we used a machine-learning clustering method to pre-determine homogenous vessels within the yellowfin tuna sub-fishery of the ETBF, expert opinion can likewise be used to identify vessels that share comparable fishing behavioural patterns. There are several important applications that will benefit from the empirical inference method we have developed. For instance, there is a need to evaluate the performance of individual fishing vessels and fleets against quantifiable targets such as bycatch performance measures or reference points, to inform management decision-making (Grafton et al., 2007; Gjertsen et al., 2010; Kirby and Ward, 2014). Our standardized bycatch rate can be used as a key indicator to measure the performance of an individual vessel/fleet relative to quantifiable targets (while also accounting for uncertainty) to identify outperforming and underperforming vessels for further investigation or corrective action. In our case study, it has allowed fishery managers to compare seabird bycatch rates of individual vessels and the fleet relative to the Australian TAP maximum permissible bycatch rate of 0.05 birds per 1000 hooks and quantitatively measure how individual vessels are performing relative to the fleet average. This can also be updated regularly to ensure responsiveness to changes in the status of bycatch species or reference points. Our inference method also allows a hierarchy of the homogenous fleet to be developed in a risk management context to prioritize resourcing and inform management decision-making. Decision rules can then be formulated based on each level of the hierarchy if considered prudent. We define three hierarchical levels based on the standardized bycatch rates (i.e. risk to seabirds), uncertainty and pre-existing management objectives (e.g. TAP: 0.05 seabirds per 1000 hooks). The “low-risk element” (i.e. those vessels with standardized bycatch rates and confidence intervals below the pre-existing limit reference point) would be considered best practice in the fishery and outperforming vessels, from which further information could be sought to determine their success in deploying mitigation measures and reducing bycatch. The “high-risk element” (i.e. those vessels with standardized bycatch rates and confidence intervals above the pre-existing limit reference point) would be considered poor-performing and prioritized for the investigation to determine what corrective action or mitigation measures are required to improve performance. The “uncertain risk element” (i.e. those vessels standardized bycatch rates above or below the pre-existing limit reference point but with confidence intervals that encompass the pre-existing limit reference point) is prioritized for further analysis to identify if their fishing operations share practices that reflect vessels in the “high-risk element”. If similar practices are identified, corrective actions can be implemented. If the analysis remains inconclusive, these vessels may be prioritized for more intensive monitoring to rapidly acquire informative data before any decision could be made about their performance. In the absence of a pre-defined bycatch performance measure, the standardized bycatch rate of the fleet could contribute to the formation of an appropriate performance measure (e.g. limit reference point) for an individual bycatch species. Conventionally, a limit reference point is defined as the level at which the risk of recruitment impairment is regarded as unacceptably high, or the minimum acceptable level of bycatch at which the measures being adopted are likely to be having the desired conservation effect (Tuck, 2011; Moore et al., 2013; DAWR, 2018). When set as a performance measure (e.g. the Australian TAP for seabirds), it provides guidance on expected levels of performance for industry and provides the means for decision-makers to evaluate and improve bycatch mitigation (Grafton et al., 2007). It also represents a uniform control limit for vessels that will drive adaptation and facilitate the robust assessment of mitigation technologies (Komoroske and Lewison, 2015). In the absence of information to determine population abundance using conventional assessments, this type of analysis can allow different stakeholders or interest groups to discuss appropriate limit reference points, which could be readily adjusted upon application or if new information on population abundance becomes available. Moreover, it can be applied in the context of “continuous improvement” until a limit reference point is defined with the objective of continually lowering the standardized bycatch rate of the fleet. The ability to use a standardized bycatch rate to measure annually the individual and fleet performance against the limit reference point can create incentives for industry to be more individually accountable of their bycatch. This can be achieved by decision-makers introducing penalties (and/or rewards) for vessels that exceed (or maintain their bycatch below) the limit reference point (Diamond, 2004; Pascoe et al., 2010). These market-based incentives could be in the form of restricting access to certain fishing areas, temporary loss of right of access and/or fines, creating a cost for sub-standard performance that would induce fishers to make choices that reduce bycatch (Diamond, 2004; Pascoe et al., 2010). This is not too dissimilar from the system of dolphin mortality limits established to manage dolphin bycatch in the purse-seine tuna fisheries of the eastern Pacific Ocean managed under the Agreement on the International Dolphin Conservation Programme (Anon, 1999; Gjertsen et al., 2010). Under this programme, a total annual limit of 5000 dolphins is set for the fishery in the Agreement Area and an equal share of this limit assigned to each applicable vessel (Anon, 1999). If at any time a vessel exceeds their dolphin mortality limit, they must cease fishing for tuna in association with dolphins, creating an incentive for improved bycatch mitigation. There is also a similar programme for the management of New Zealand sea lion (Phocarctos hookeri) mortalities in the New Zealand squid fishery, with a fishing-related mortality limit derived from a Bayesian model (Breen et al., 2003) set annually (Chilvers, 2008). Once the limit is reached within a season, the fishery is then closed, creating an incentive for fishers to reduce their bycatch (Robertson and Chilvers, 2011). While our standardized bycatch rate cannot be used to measure current population status (initial or current abundance), it can be used to monitor the performance of individual vessels and the fleet relative to the performance measure for an individual species. Of course, this assumes that decision-makers have access to data at a species taxonomic level that can be trusted. Fisher-reported logbook data have often been found to be inaccurate and inconsistent with at-sea observer data from the same trip, due to fishers either misreporting, under-reporting, over reporting, or non-reporting their bycatch (Sampson, 2011; Mangi et al., 2016; Macbeth et al., 2018). While in this case study we used logbook data that have been verified (using an electronic monitoring programme) (Emery et al., 2019a, b), our model is not constrained to fisheries with verifiable logbook data. It can easily be applied to fisheries with unverified logbook data or extrapolated at-sea observer data (assuming coverage is sufficient) but noting the issues and caveats with precision remain the same as if an alternative model was run using that data (Wakefield et al., 2018). We developed a model to estimate standardized individual vessel and fleet bycatch rates that can be widely applied, is simple and accessible for fisheries with limited data, can deal with uncertainty in rate estimation, and can be easily interpreted in a risk context. Risk-based approaches or frameworks are useful for decision-makers to prioritize scarce resources (both in terms of further investigation or corrective action). Our model can also be readily updated to determine whether a vessel’s bycatch rate changes over time or following intervention and has the potential to include additional information such as location and seasonality as covariates. Lastly, this approach could be tailored to each bycatch issue or situation and combined with additional risk-based models, such as fisheries compliance risk assessments (e.g. AFMA, 2017), to provide a more comprehensive risk framework for the fishery. Acknowledgements We would like to acknowledge the Australian Fisheries Management Authority (AFMA) for providing the commercial catch and effort data for the Eastern Tuna and Billfish Fishery (ETBF). We also would like to thank the AFMA ETBF manager Don Bromhead for fruitful discussions on the analysis and manuscript. Lastly, we thank Rupert Summerson (ABARES) for producing the map of the fishery. References AFMA. 2017 . National Compliance 2017-19 Risk Assessment Methodology. Australian Fisheries Management Authority, Canberra. Anon. 1999 . Agreement on the International Dolphin Conservation Programme (Amended). p. 23 . Basu A. P. , Rigdon S. E. 1986 . Examples of parametric empirical Bayes methods for the estimation of failure processes for repairable systems. In Reliability and Quality Control , pp. 47 – 55 . Ed. by Basu A. P. Elsevier, Amsterdam . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Benoît H. P. , Allard J. 2009 . Can the data from at-sea observer surveys be used to make general inferences about catch composition and discards? Canadian Journal of Fisheries and Aquatic Sciences , 66 : 2025 – 2039 . Google Scholar Crossref Search ADS WorldCat Bishop J. , Venables W. N. , Dichmont C. M. , Sterling D. J. 2008 . Standardizing catch rates: is logbook information by itself enough? ICES Journal of Marine Science , 65 : 255 – 266 . Google Scholar Crossref Search ADS WorldCat Boonstra W. J. , Hentati-Sundberg J. 2016 . Classifying fishers’ behaviour. An invitation to fishing styles . Fish and Fisheries , 17 : 78 – 100 . Google Scholar Crossref Search ADS WorldCat Breen P. A. , Hilborn R. , Maunder M. N. , Kim S. W. 2003 . Effects of alternative control rules on the conflict between a fishery and a threatened sea lion (Phocarctos hookeri) . Canadian Journal of Fisheries and Aquatic Sciences , 60 : 527 – 541 . Google Scholar Crossref Search ADS WorldCat Carlin, B. P. and Louis, T. A. 2009. Bayesian Methods for Data Analysis. CRC Press Chapman and Hall, Florida, USA. Carlin B. P. , Gelfand A. E. 1991 . A sample reuse method for accurate parametric empirical bayes confidence intervals . Journal of the Royal Statistical Society. Series B (Methodological) , 53 : 189 – 200 . Google Scholar Crossref Search ADS WorldCat Chilvers B. L. 2008 . New Zealand sea lions Phocarctos hookeri and squid trawl fisheries: bycatch problems and management options . Endangered Species Research , 5 : 193 – 204 . Google Scholar Crossref Search ADS WorldCat Commonwealth of Australia. 2018 . Threat Abatment Plan for the incidental catch (or bycatch) of seabirds during oceanic longline operations. Department of Environment and Energy, Canberra. https://www.legislation.gov.au/Details/F2018L01562 (last accessed 2 November 2019). Cortés V. , Arcos J. M. , González-Solís J. 2017 . Seabirds and demersal longliners in the northwestern Mediterranean: factors driving their interactions and bycatch rates . Marine Ecology Progress Series , 565 : 1 – 16 . Google Scholar Crossref Search ADS WorldCat DAWR. 2018 . Guidelines for the Implementation of the Commonwealth Fisheries Bycatch Policy. Department of Agriculture and Water Resources , Canberra . 52 .pp. OpenURL Placeholder Text WorldCat Diamond S. L. 2004 . Bycatch quotas in the Gulf of Mexico shrimp trawl fishery: can they work? Reviews in Fish Biology and Fisheries , 14 : 207 – 237 . Google Scholar Crossref Search ADS WorldCat Duarte R. , Azevedo M. , Afonso-Dias M. 2009 . Segmentation and fishery characteristics of the mixed-species multi-gear Portuguese fleet . ICES Journal of Marine Science , 66 : 594 – 606 . Google Scholar Crossref Search ADS WorldCat Emery T. J. , Noriega R. , Williams A. J. , Larcombe J. 2019a . Changes in logbook reporting by commercial fishers following the implementation of electronic monitoring in Australian Commonwealth fisheries . Marine Policy , 104 : 135 – 145 . Google Scholar Crossref Search ADS WorldCat Emery T. J. , Noriega R. , Williams A. J. , Larcombe J. 2019b . Measuring congruence between electronic monitoring and logbook data in Australian Commonwealth longline and gillnet fisheries . Ocean & Coastal Management , 168 : 307 – 321 . Google Scholar Crossref Search ADS WorldCat Faunce C. H. , Barbeaux S. J. 2011 . The frequency and quantity of Alaskan groundfish catcher-vessel landings made with and without an observer . ICES Journal of Marine Science , 68 : 1757 – 1763 . Google Scholar Crossref Search ADS WorldCat Gardner B. , Sullivan P. J. , Epperly S. , Morreale S. J. 2008 . Hierarchical modeling of bycatch rates of sea turtles in the western North Atlantic . Endangered Species Research , 5 : 279 – 289 . Google Scholar Crossref Search ADS WorldCat Gilman E. , Clarke S. , Brothers N. , Alfaro-Shigueto J. , Mandelman J. , Mangel J. , Petersen S. et al. . 2008 . Shark interactions in pelagic longline fisheries . Marine Policy , 32 : 1 – 18 . Google Scholar Crossref Search ADS WorldCat Gjertsen H. , Hall M. , Squires D. 2010 . Incentives to address bycatch issues. In Conservation and Management of Transnational Tuna Fisheries , pp. 225 – 248 . Ed. by Allen R. , Joseph J. , Squires D. . Blackwell Publishing, Iowa, USA . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Grafton Q. , Kompas T. , McLoughlin R. , Rayns N. 2007 . Benchmarking for fisheries governance . Marine Policy , 31 : 470 – 479 . Google Scholar Crossref Search ADS WorldCat Greenwood M. , Yule G. U. 1920 . An inquiry into the nature of frequency distributions representative of multiple happenings with particular reference to the occurrence of multiple attacks of disease or of repeated accidents . Journal of the Royal Statistical Society , 83 : 255 – 279 . Google Scholar Crossref Search ADS WorldCat Hatch J. M. 2018 . Comprehensive estimates of seabird–fishery interactions for the US Northeast and mid-Atlantic . Aquatic Conservation: Marine and Freshwater Ecosystems , 28 : 182 – 193 . Google Scholar Crossref Search ADS WorldCat Holley J.-F. , Marchal P. 2004 . Fishing strategy development under changing conditions: examples from the French offshore fleet fishing in the North Atlantic . ICES Journal of Marine Science , 61 : 1410 – 1431 . Google Scholar Crossref Search ADS WorldCat Khakzad N. , Khan F. , Paltrinieri N. 2014 . On the application of near accident data to risk analysis of major accidents . Reliability Engineering & System Safety , 126 : 116 – 125 . Google Scholar Crossref Search ADS WorldCat Kirby D. S. , Ward P. 2014 . Standards for the effective management of fisheries bycatch . Marine Policy , 44 : 419 – 426 . Google Scholar Crossref Search ADS WorldCat Komoroske L. M. , Lewison R. L. 2015 . Addressing fisheries bycatch in a changing world . Frontiers in Marine Science , 2 : 1 – 11 . Google Scholar Crossref Search ADS WorldCat Larcombe J. , Noriega R. , Timmiss T. 2016 . Catch reporting under e-monitoring in the Australian Pacific longline fishery. Second Reporting and E-Monitoring Intersessional Working Group Meeting. WCPFC-2016 ERandEMWG2-DPO1, 20 pp. https://www.wcpfc.int/node/27541 (last accessed 2 November 2019). Lent R. , Squires D. 2017 . Reducing marine mammal bycatch in global fisheries: an economics approach . Deep Sea Research Part II: Topical Studies in Oceanography , 140 : 268 – 277 . Google Scholar Crossref Search ADS WorldCat Macbeth W. G. , Butcher P. A. , Collins D. , McGrath S. P. , Provost S. C. , Bowling A. C. , Geraghty P. T. et al. . 2018 . Improving reliability of species identification and logbook catch reporting by commercial fishers in an Australian demersal shark longline fishery . Fisheries Management and Ecology , 25 : 186 – 202 . Google Scholar Crossref Search ADS WorldCat Mangi S. C. , Smith S. , Catchpole T. L. 2016 . Assessing the capability and willingness of skippers towards fishing industry-led data collection . Ocean & Coastal Management , 134 : 11 – 19 . Google Scholar Crossref Search ADS WorldCat Martin S. L. , Stohs S. M. , Moore J. E. 2015 . Bayesian inference and assessment for rare-event bycatch in marine fisheries: a drift gillnet fishery case study . Ecological Applications , 25 : 416 – 429 . Google Scholar Crossref Search ADS PubMed WorldCat Moore J. E. , Curtis K. A. , Lewison R. L. , Dillingham P. W. , Cope J. M. , Fordham S. V. , Heppell S. S. et al. . 2013 . Evaluating sustainability of fisheries bycatch mortality for marine megafauna: a review of conservation reference points for data-limited populations . Environmental Conservation , 40 : 329 – 344 . Google Scholar Crossref Search ADS WorldCat Myers R. A. , Barrowman N. J. , Hilborn R. , Kehler D. G. 2002 . Inferring Bayesian priors with limited direct data: applications to risk analysis . North American Journal of Fisheries Management , 22 : 351 – 364 . Google Scholar Crossref Search ADS WorldCat Pascoe S. , Innes J. , Holland D. , Fina M. , Thébaud O. , Townsend R. , Sanchirico J. et al. . 2010 . Use of incentive-based management systems to limit bycatch and discarding . International Review of Environmental and Resource Economics , 4 : 123 – 161 . Google Scholar Crossref Search ADS WorldCat Patterson H. , Larcombe J. , Nicol S. , Curtotti R. 2018 . Fishery Status Reports 2018. Australian Bureau of Agricultural and Resource Economics and Sciences, Canberra. Pelletier D. , Ferraris J. 2000 . A multivariate approach for defining fishing tactics from commercial catch and effort data . Canadian Journal of Fisheries and Aquatic Sciences , 57 : 51 – 65 . Google Scholar Crossref Search ADS WorldCat Quigley J. , Hardman G. , Bedford T. , Walls L. 2011 . Merging expert and empirical data for rare event frequency estimation: pool homogenisation for empirical Bayes models . Reliability Engineering & System Safety , 96 : 687 – 695 . Google Scholar Crossref Search ADS WorldCat Robertson B. C. , Chilvers B. L. 2011 . The population decline of the New Zealand sea lion Phocarctos hookeri: a review of possible causes . Mammal Review , 41 : 253 – 275 . Google Scholar Crossref Search ADS WorldCat Sampson D. B. 2011 . The accuracy of self-reported fisheries data: Oregon trawl logbook fishing locations and retained catches . Fisheries Research , 112 : 59 – 76 . Google Scholar Crossref Search ADS WorldCat Suuronen P. , Gilman E. 2019 . Monitoring and managing fisheries discards: new technologies and approaches . Marine Policy , 103554 . OpenURL Placeholder Text WorldCat Tuck G. N. 2011 . Are bycatch rates sufficient as the principal fishery performance measure and method of assessment for seabirds? Aquatic Conservation: Marine and Freshwater Ecosystems , 21 : 412 – 422 . Google Scholar Crossref Search ADS WorldCat Wakefield C. B. , Hesp S. A. , Blight S. , Molony B. W. , Newman S. J. , Hall N. G. 2018 . Uncertainty associated with total bycatch estimates for rarely-encountered species varies substantially with observer coverage levels: informing minimum requirements for statutory logbook validation . Marine Policy , 95 : 273 – 282 . Google Scholar Crossref Search ADS WorldCat Walsh W. A. , Ito R. Y. , Kawamoto K. E. , McCracken M. 2005 . Analysis of logbook accuracy for blue marlin (Makaira nigricans) in the Hawaii-based longline fishery with a generalized additive model and commercial sales data . Fisheries Research , 75 : 175 – 192 . Google Scholar Crossref Search ADS WorldCat Walsh W. A. , Kleiber P. , McCracken M. 2002 . Comparison of logbook reports of incidental blue shark catch rates by Hawaii-based longline vessels to fishery observer data by application of a generalized additive model . Fisheries Research , 58 : 79 – 94 . Google Scholar Crossref Search ADS WorldCat © International Council for the Exploration of the Sea 2020. All rights reserved. For permissions, please email: [email protected] This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)