Enhancing the utility of known-biomass production models: a case study of the Bay of Biscay and Iberian Coast ecoregion

Case study: modelling commercial stocks in the Bay of Biscay and Iberian Coast ecoregion

The study area includes the ICES subareas 3-4 and 6-9, covering the Iberian Peninsula and Bay of Biscay ecoregion and the Celtic Sea ecoregion (ICES 2021). Overall, these ecoregions are characterized by marked seasonal mixing and a stratification of water masses typical of temperate seas. Modifications to this pattern may be due to wind-driven upwelling, river runoff and tidal processes, which increase the system’s productivity with large variations across the region (ICES 2021ICES 2021. Bay of Biscay and the Iberian Coast ecoregion – Ecosystem Overview. ICES Advice: Ecosystem Overviews. Report. https://doi.org/10.17895/ices.advice.9436 ). Habitats far from the coast are shaped by the influence of Atlantic waters in the Bay of Biscay and western Iberia (Abad et al. 2020Abad E., Pennino M.G., Valeiras J., et al. 2020. Integrating spatial management measures into fisheries: The Lepidorhombus spp. case study. Mar. Policy 116: 103739. https://doi.org/10.1016/j.marpol.2019.103739 ).

The target species studied were megrim (Lepidorhombus whiffiagonis), white anglerfish (Lophius piscatorius) and European hake (Merluccius merluccius). These species are generally targeted together in a mixed fishery and have traditionally been considered of great commercial importance for fishing because of their high economic value.

All species are widely distributed throughout Atlantic shelf waters, and the selected stocks are assessed by the ICES in two units, northern and southern stocks, which are the ones we used in our study.

Dataset

Data were obtained from the ICES database (https://www.ices.dk) through the icesSAG package (Colin et al. 2022Colin M., Scott L., Arni M., Pinto C. 2022. icesSAG: Stock Assessment Graphs Database Web Services. R package version 1.4.0. https://CRAN.R-project.org/package=icesSAG ) of R software (R Core Team 2021R Core Team 2021. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/ ). Indeed, this package allows access to the web services of the ICES database and download of the required data directly in R. For this study, spawning stock biomass (SSB, in tonnes) and annual catches (in tonnes) were collected for each stock for the time period reported in Table 1.

In order to test the potential effect of environmental variability on stock dynamics, we used two indices: the Atlantic Multidecadal Oscillation (AMO) and the North Atlantic Oscillation (NAO). The AMO is a long-term climatic variability index, which is defined as the North Atlantic sea surface temperature anomaly between Greenland and the Equator. The AMO provides an ongoing series of long-duration changes in sea surface temperature of the North Atlantic Ocean, with cool and warm phases that may last for 20 to 40 years at a time (Kerr 2000Kerr R.A. 2000. A North Atlantic climate pacemaker for the centuries. Science 288 (5473): 1984-1986. https://doi.org/10.1126/science.288.5473.1984 , Enfield et al. 2001Enfield D., Mestas-Nunez A., Trimble P. 2001. The Atlantic multidecadal oscillation and its relation to rainfall and river flows in the continental U.S. Geophys. Res. Lett. 28: 2077-2080. https://doi.org/10.1029/2000GL012745 ). Its inclusion allows us to address longer-term variations that can influence fish population dynamics. We also included the NAO index in our study because it is one of the major sources of climate variability in the northern hemisphere (Hurrell 1995Hurrell J.W. 1995. Decadal trends in the North Atlantic Oscillation: Regional Temperatures and Precipitation. Science 269: 676-679. https://doi.org/10.1126/science.269.5224.676 ). The NAO serves as a comprehensive representation of the winter maritime climate (e.g. swell) and the trajectories of storms as they pass over the North Atlantic and therefore has a major influence on wind and precipitation in our region (Hurrell and Deser 2009Hurrell J.W., Deser C. 2009. North Atlantic climate variability: the role of the North Atlantic Oscillation. J. Mar. Syst. 78: 28-41. https://doi.org/10.1016/j.jmarsys.2008.11.026 ). We obtained the NAO and AMO annual averages from the NOAA (available at https://www.cpc.ncep.noaa.gov/data/teledoc/telecontents.shtml).

Data analyses

For a better understanding of the stock dynamics of the three demersal species considered in the study, the proposed analytical strategy involved four steps: (1) annual SP was calculated and KBPMs were fitted to each of the six stock units listed in Table 1; (2) for the stocks that showed possible production changes in the individual KBPM fit, a retrospective analysis was performed to determine whether such changes had occurred for reasons other than biomass variability; (3) multispecies KBPMs were implemented to analyse the dynamics of the southern and northern communities of the Bay of Biscay and Iberian Coast ecoregion, and also of the global component aggregating both areas; and (4) the AMO and NAO indices were included in KBPMs to assess the environmental influence on stock status for the stocks identified in (2) as experiencing production change caused by factors other than biomass variability.

Known-biomass production models

The only data requirements for SP models with known biomass are the time series for catch and biomass (estimated values of biomass produced by a data-rich stock assessment model). Traditional SP models estimate biomass as “exploitable biomass”, that is, as the fraction of biomass accessible to fishing gears based on their selectivity. For production models with a known biomass, estimates of exploitable biomass are not available, but total biomass is available in some cases and SSB is always available. Hence, annual SP is calculated for each stock as:

for t=1, .., T, where T is the last year of the catch time series, $SS{B}_t$ the spawning biomass at the beginning of year t, and $C_t$ the catch during year t. Then, the average annual spawning biomass $\overline{SS{B}_t}$ that produced $S{P}_t$ is computed as $\overline{SS{B}_t}$ = (SS $B_t$ +SS $B_{t+1}$ )/2. The surplus production-biomass relationship is fitted based on $\overline{SS{B}_t}$ and $S{P}_t$ as follows:

for t=1, .., T, where a, b, and $\nu$ are model parameters and ${\epsilon }_t$ are model residuals that are assumed to be normally distributed. The parameter $\nu$ is the shape parameter introduced by Pella and Tomlinson (1969)Pella J. J., Tomlinson P. K. 1969. A generalized stock-production model. Bull I-ATTC 13: 421-58. to allow asymmetry in the SP curve. The asymmetry parameter, $\nu$ , is difficult to estimate accurately, so henceforth we assume $\nu$ =2, which corresponds to the symmetric SP curve of the Schaefer formulation (1957)Schaefer M. B. 1957. A study of the dynamics of the fishery for yellowfin tuna in the Eastern Tropical Pacific Ocean. Bull I-ATTC 2: 247-285..

The BRP estimates derived from model (2) are $SS{B}_0$ =-a/b (where $SS{B}_0$ is the maximum value of SSB when SP=0, which is the maximum unfished equilibrium SSB), SS $B_{msy}$ = $SS{B}_0$ /2 (where $SS{B}_{msy}$ is the value of SSB when SP=MSY), MSY =aSS $B_{msy}$ +b $SS{B}_{msy}^2$ and $F_{msy}$ = a/2. The parameter a is algebraically equivalent to the intrinsic population growth rate in the logistic population growth model (Prager 1992Prager M.H. 1992. ASPIC: A Surplus-Production Model Incorporating Covariates. Coll. Vol. Sci. Pap. Int. Comm. Conserv. Atl. Tunas (ICCAT) 28: 218-229., Jacobson et al. 2002Jacobson L., Cadrin S., Weinberg J. 2002. Tools for Estimating Surplus Production and FMSY in Any Stock Assessment Model. N. Am. J. Fish. Manag. 22: 326-338. https://doi.org/10.1577/1548-8675(2002)022<0326:TFESPA>2.0.CO;2 ). Finally, the data-rich estimates of SSB and the exploitation rate in the final year were compared with the reference points derived from KBPMs to draw conclusions regarding stock status. These results were crosschecked with the conclusions obtained from the data-rich stock BRPs.

Retrospective analysis

In order to identify possible changes in SP caused by reasons other than biomass variability, a retrospective analysis was performed for the stocks that showed possible population changes in the individual KBPM fits. Retrospective patterns are usually defined as systematic changes in the estimates of assessment model-derived quantities that occur as additional years of data are removed from a stock assessment (Hurtado-Ferro et al. 2015Hurtado-Ferro F., Szuwalski C.S., Valero J.L., et al. 2015. Looking in the rear-view mirror: bias and retrospective patterns in integrated, age-structured stock assessment models. ICES Mar. Sci. 72: 99-110. https://doi.org/10.1093/icesjms/fsu198 ). Here, we performed a retrospective analysis by removing three years at a time (2020-2017-2014-2011-2008-2005) for northern hake (northern stock 7-8abd), and assessing whether the model fit changed in relation to parameters a and b. We focus on this stock because, as we will see in the KBPM fit results, it is the one which exhibits a significant, abrupt and sustained change in SP and SSB.

Multispecies KBPMs

Ecosystem-based fisheries management (EBFM) has been proposed as a place-based rather than a species-based management approach (Fogarty 2014Fogarty M.J. 2014. The art of ecosystem-based fishery management. Can. J. Fish. Aquat. Sci. 71: 479-490. https://doi.org/10.1139/cjfas-2013-0203 ) that can account for fishery effects on food web structure and function, biodiversity balance and yield objectives, in addition to improving management in systems characterized by large seasonal, decadal and long-term climate variability, such as the Atlantic Ocean.

One possibility of KBPMs is that they can be easily fitted as multispecies stock assessment models by aggregating each stock biomass and catch time series per year (Bundy et al. 2012Bundy A., Bohaboy E.C., Hjermann D.O., et al. 2012. Common patterns, common drivers: comparative analysis of aggregate surplus production across ecosystems. Mar. Ecol. Prog. Ser. 459: 203-218. https://doi.org/10.3354/meps09787 ). This possibility could allow assessment of the species community, defining ecosystem-level reference points and identifying the ecosystem status following EBFM requirements (Mueter and Megrey 2006Mueter F., Megrey B. 2006. Using multi-species surplus production models to estimate ecosystem-level maximum sustainable yields. Fish. Res. 81: 189-201. https://doi.org/10.1016/j.fishres.2006.07.010 ). Here, we performed a multispecies approximation, aggregating southern stocks (in Iberian Peninsula Atlantic waters) and northern stocks (northwards) into two components so as to assess these two parts of the ecosystems separately. Additionally, northern and southern stocks were also aggregated into a single component to assess the status of the demersal community.

Despite its simplicity, the KBPM multispecies approach provides estimates of multispecies MSY reference points that are more appropriate than the sum of single-species MSYs, providing an alternative approach to setting limits on overall removals from an exploited species complex (Mueter and Megrey 2006Mueter F., Megrey B. 2006. Using multi-species surplus production models to estimate ecosystem-level maximum sustainable yields. Fish. Res. 81: 189-201. https://doi.org/10.1016/j.fishres.2006.07.010 ).

Testing environmental variability on stock status

In addition to modelling SP as a function of spawning biomass, we also investigated the relationships between SP and the NAO and AMO environmental indices in order to test whether productivity changes in response to environmental fluctuations for the stocks identified as undergoing a production change caused by reasons other than biomass variability. Environmental effects are modelled as either additive effects (Bundy et al. 2012Bundy A., Bohaboy E.C., Hjermann D.O., et al. 2012. Common patterns, common drivers: comparative analysis of aggregate surplus production across ecosystems. Mar. Ecol. Prog. Ser. 459: 203-218. https://doi.org/10.3354/meps09787 ) or multiplicative effects (Mueter and Megrey 2006Mueter F., Megrey B. 2006. Using multi-species surplus production models to estimate ecosystem-level maximum sustainable yields. Fish. Res. 81: 189-201. https://doi.org/10.1016/j.fishres.2006.07.010 ) as follows:

for t=1, .., T, where $X_{t-lag}$ is the environmental covariate measured in year $t-lag$ , where $lag$ is a value selected from the sequence 1, …, $n_{lag}$ according to the highest Pearson correlation among the $S{P}_t$ and $X_{t-lag}$ time series, and $n_{lag}$ is the maximum number of lags for which it is considered that a correlation may exist. More precisely, Pearson’s correlation between $X_{t-lag}$ and $S{P}_t$ time series for $lag$ from 1 to $n_{lag}$ =4 was computed, and then the $lag$ was finally set equal to the one corresponding to the highest correlation.

The models were fitted using R software (R Core Team 2021R Core Team 2021. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/ ). To fit the non-linear models, we used the nls (non-linear least square) function. The R code for the analytical strategy is available at this GitHub repository https://github.com/FranIzquierdo (https://github.com/IMPRESSPROJECT/Known-biomass-production-models-the-case-study-of-the-Bay-the-Biscay-and-Iberian-Coast-Ecoregion).

Predicted values derived from the original fit without environmental covariables, Equation (2), and from both environmental fits were compared with the observed SP values through the following error measures: mean absolute percentage error (MAPE) and root mean square error (RMSE).

where T is the total length of the time SP series, ${\widehat{ SP}}_t$ the predicted SP at time t and $\left|•\right|$ the absolute value function. In addition, we also compared performance of the models according to the Akaike information criterion (AIC).

where k is the number of estimated parameters in the model and $\widehat{L}$ the maximized value of the likelihood function for the model.

Known-biomass production models

The average annual spawning biomass, annual SP and model fit for these data were computed for each individual stock (Fig. 1). According to Walters et al. (2008)Walters C.J., Hilborn R., Christensen V. 2008. Surplus production dynamics in declining and recovering fish populations. Can. J. Fish. Aquat. Sci. 65: 2536-2551. https://doi.org/10.1139/F08-170 , several patterns can be identified based on a visual inspection of $\overline{SS{B}_t}$ and $S{P}_t$ pairs: (i) stationary (production varying randomly around a mean production curve); (ii) upward or clockwise hooks-cycles (increasing production and stock followed by a decrease in production and stock, or decreasing production and stock followed by an increase in stock and production); and (iii) downward or counterclockwise hooks-cycles (higher production rates during stock declines than during subsequent recoveries).

Northern anglerfish and both hake stocks (southern and northern) show clockwise hooks-cycles. More precisely, when biomass decreases, SP is below expected; then, when SP surpasses what is expected, biomass increases again. A counterclockwise hooks-cycle is observed for southern anglerfish; higher production rates during SSB decline were detected from 1980 to 1994, then the general pattern was an increase in SSB and a decrease in SP.

The $\overline{SS{B}_t}$ and $S{P}_t$ pairs of the northern and southern megrim correspond to the left side of the symmetric Schaefer curve, which leads us to conclude that the maximum SP in these stocks is estimated with high uncertainty. A clockwise hooks-cycle was observed for northern and southern megrim: when SP was below expected, the SSB decreased until SP crossed the expected value and then SSB increased.

The KBPM fit information (Fig. 1) is completed by the biological reference point estimates derived in Table 2 and the parameter estimates included in the Supplementary Material (Table S1). The comparison of the last year of SSB and F (exploitation rate derived from the prior data-rich stock assessment fit) with the corresponding KBPM reference points (SSB_msy and F_msy) allows us to discuss the actual stock status compared with BRP independently of a specific stock recruitment function. Additionally, the conclusions derived from the KBPM BRPs can be crosschecked with the data-rich stock assessment outcomes. As mentioned above, this information can be easily obtained using the icesSAG package (code available on GitHub repository). For southern and northern megrims and anglerfish, the F value in the last year was below F_msy for both estimates (data-rich and KBPM), leading to the same conclusion about their exploitation status. However, although both lead to the same final conclusions, it is worth mentioning that for anglerfish the KBPM’s F_msy is considerably higher than the data-rich one, whereas for megrims both estimates were closer to each other. For northern hake the conclusions obtained from the two approaches do not coincide because the data-rich model concluded that stock is exploited at F_msy, whereas KBPM concluded that fishing pressure is considerably below F_msy. Finally, for the southern hake stock the F value in the last year was above the data-rich F_msy but not above the KBPM F_msy value, leading to different conclusions about the fishing pressure applied to this stock.

The MSY/SSB₀ ratio is informative about current stock productivity, and according to the values in Table 2, megrim stocks can be classified as the least productive stocks of those studied, whereas hake (northern and southern) and northern anglerfish are the most productive.

The substantial, abrupt and persistent change in SP and SSB detected in northern hake stock can lead to the hypothesis that, in addition to biomass variability, there are other reasons for this shift, which is why we performed a retrospective analysis of this stock.

Retrospective analysis

A retrospective analysis was performed by removing three years at a time (2020-2017-2014-2011-2008-2005) for northern hake stock (northern stock 7-8abd). Figure 2 reports on the retrospective analysis, Table 3 provides the corresponding reference points, and the Supplementary Material (Table S2) provides the parameter estimates.

The first and second retrospective fits (end of years 2017 and 2014) showed a slight increase in MSY, which was clearer in the third retrospective fit (end of year 2011), with a value of 136261 (an increase of 16% in relation to the value in the 1978-2020 fit), whereas stock productivity remained stable (MSY/SSB₀ around 0.3) and SSB₀ and SSB_msy increased. Conversely, the fourth and fifth retrospective fits (end of years 2008 and 2005) showed an increase in productivity (around 60% in relation to the original MSY/SSB₀ value), whereas the corresponding SSB₀ and SSB_msy decreased by around 68% in relation to the original fit values (end of year 2020). The retrospective fits lead to the conclusion that from 2005 to 2017 the carrying capacity increased, whereas the relative productivity decreased even though the SSB increased. Consequently, this retrospective analysis of the production curve using KBPMs identifies a pattern in the dynamics of the northern hake stock that requires thorough examination to uncover its underlying causes. Step (4) of our analytical approach aimed to investigate whether these variations could be attributable to environmental fluctuations.

Multispecies KBPMs

We conducted a multispecies approximation, aggregating the northern and southern stocks into two components to be assessed separately (Fig. 3, Table 4 and Table S3 in the Supplementary Material). Furthermore, the global area was also analysed through a single component aggregating northern and southern stocks (Fig. 3, Table 4 and Table S3 in the Supplementary Material). This practical usage of KBPM allows the demersal community status to be analysed from an ecosystem point of view. Northern and southern demersal community patterns are clockwise hooks-cycles (decreasing production and SSB followed by an increase in both production and SSB). The global community pattern is clearly determined by the northern community because of its high size in relation to the southern one. Additionally, this simple empirical approach provides multispecies MSY BRPs (Table 4). In particular, based on the MSY/SSB₀ values, it can be concluded that the southern demersal community exhibits greater productivity than its northern counterpart (MSY/SSB₀ values of 0.32 and 0.23, respectively), with the MSY northern community value being approximately eight times that of the southern one (Table 4).

Testing the environmental variability in stock status

The relationship between SP and the NAO and AMO environmental variables was investigated for northern hake because drivers of production changes in marine fish communities can include climate forcing. The scatter plot between the values of SP and the environmental covariates at different year lags from 1 to 4 was derived and the corresponding correlation was calculated (Fig. 4). The maximum correlation value corresponded to lag 4 of the AMO covariable, so additive model and multiplicative environmental models, Equations (3) and (4), were fitted including this covariable. A significant environmental effect was found in the multiplicative model, where c=6.358e-01, whereas the environmental effect in the additive model was not significant (Table S4 in the Supplementary Material).

The predicted values derived from the original fit without the environmental covariable, Equation (2), and from both environmental fits were reported with the observed SP values (Fig. 5). Graphically, we found no clear differences, so the MAPE and RMSE measures were computed and reported for each of the models in Table 5 with the AIC values. The lowest values of AIC, MASE and MAPE corresponded to the multiplicative model; however, the differences between the error measure values of the different models were minimal. This finding leads to the conclusion that although the multiplicative model detected an environmental effect on stock production, this does not seem to be the main factor explaining the change in stock SP, or that more local environmental covariables would have to be included in the model to identify the environmental factor as a main one.