Abstract

Environmental conditions can create spatial and temporal variability in growth and distribution processes, yet contemporary stock assessment methods often do not explicitly address the consequences of these patterns. For example, stock assessments often assume that body weight-at-age (i.e. size) is constant across the stocks’ range, and may thereby miss important spatio-temporal patterns. This is becoming increasingly relevant given climate-driven distributional shifts, because samples for estimating size-at-age can be spatially unbalanced and lead to biases when extrapolating into unsampled areas. Here, we jointly analysed data on the local abundance and size of walleye pollock (Gadus chalcogrammus) in the Bering Sea, to demonstrate a tractable first step in expanding spatially unbalanced size-at-age samples, while incorporating fine-scale spatial and temporal variation for inclusion in stock assessments. The data come from NOAA’s bottom trawl survey data and were evaluated using a multivariate spatio-temporal statistical model. We found extensive variation in size-at-age at fine spatial scales, though specific patterns differed between age classes. In addition to persistent spatial patterns, we also documented year-to-year differences in the spatial patterning of size-at-age. Intra-annual variation in the population-level size-at-age (used to generate the size-at-age matrix in the stock assessment) was largely driven by localized changes in fish size, while shifts in species distribution had a smaller effect. The spatio-temporal size-at-age matrix led to marginal improvement in the stock assessment fit to the survey biomass index. Results from our case study suggest that accounting for spatially unbalanced sampling improved stock assessment consistency. Additionally, it improved our understanding on the dynamics of how local and population-level demographic processes interact. As climate change affects fish distribution and growth, integrating spatiotemporally explicit size-at-age processes with anticipated environmental conditions may improve stock-assessment forecasts used to set annual harvest limits.

Introduction

Fisheries managers depend on stock assessments to make informed decisions about conservation and management measures (Methot, 2009). State-of-the-art methods typically involve fitting an age-structured population dynamics model to evaluate the current and projected status of the resource including biomass and fishing mortality (Hilborn and Walters, 1992). The population dynamics model relies on data from surveys and fisheries to estimate survival patterns and age-specific demographic rates, including growth, mortality, recruitment, and maturation. Key stock characteristics estimated from assessment models are sensitive to changes in these demographic rates (Thorson et al., 2015a). Therefore, the accuracy and precision of the data ultimately impact the value of stock assessments for managers, especially when assumptions such as constant growth rates can be avoided (Methot, 2009; Francis, 2011).

Survey and fishery data on catch rates and composition inevitably have biases and gaps that present challenges to using them in stock assessments (Thorson, 2014; Breivik et al., 2021; Ducharme-Barth et al., 2022). While well-designed fishery-independent surveys attempt to account for known sources of variation with standardized protocols (Gunderson, 1993), there are inevitable disruptions in survey implementation and unobservable sources of variation that make this challenging (Kimura and Somerton, 2007). For instance, sampling intensity may not be equally distributed across a species’ range (Thorson, 2014), there may be aberrant extreme catches (Thorson et al., 2011), and the catchability of survey gear may vary (Kotwicki and Ono, 2018). Resource limitations and weather conditions may contribute to unforeseen gaps in spatial coverage (Adams et al., 2021). In 2020, for instance, many planned surveys were canceled due to COVID-19 restrictions (e.g. NOAA, 2020a,b, c, d). As fish distributions shift due to climate change (Mueter et al., 2011; Pinsky et al., 2020, Rooper et al., 2021), available data may increasingly fail to align with species’ actual ranges (Link et al., 2010; O'Leary et al., 2020) and will exacerbate the challenge of handling biases and gaps in data.

Spatio-temporal models are being used more commonly to address these issues (Maunder and Punt, 2004; Francis, 2011; Maunder et al., 2020; O'Leary et al., 2020; Breivik et al., 2021; Ducharme-Barth et al., 2022). By accounting for persistent spatial differences (i.e. spatial variation) and spatial differences that change over time (i.e. spatio-temporal variation) in the catch rate data, these models can provide better estimates of annual abundance in areas with missing data (e.g. Breivik et al., 2021), extrapolate to unsampled years (e.g. O'Leary et al., 2020), and correct for differences in sampling intensity across the range (e.g. Maunder et al., 2020). Spatio-temporal models have begun to similarly be applied to age composition data (Thorson, 2014; Thorson and Haltuch, 2018; O'Leary et al., 2020). These models provide a way to predict fish density by standardizing catch rates by the area sampled (“area-weighting”) and to predict composition information, such as age and size, by standardizing composition data by the fish density (“density-weighting”) of each location (Thorson et al., 2020). Spatio-temporal models were found to improve prediction and accuracy of abundance estimates over other standardization methods (Thorson et al.,2015b; Gruss et al., 2019; Zhou et al., 2019).

However, spatio-temporal models have rarely been applied to standardizing size-at-age (here defined as weight-at-age), another critical component of assessment models. While there are different methods used to incorporate growth rates in stock assessment models (Helser and Brodziak, 1998; Clark and Hare, 2002; Minte-Vera, 2004; Whitten et al., 2013; Thorson and Minte-Vera, 2016), often an annual size-at-age is derived from size measurements of survey or fishery data in an empirical size-at-age matrix (e.g. Kuriyama et al., 2016; Ianelli et al., 2021). Standardizing data for size-at-age is inherently more complicated than standardizing catch rate for abundance indices because size data results from the interplay between variation in two processes: growth (and resulting size) within different habitats, and the stock abundance in those habitats. First, environmental differences in food availability, temperature, or oxygen may impact local growth rates, resulting in variation in local size-at-age (DeVries and Frie, 1996). Second, fish populations are distributed in space via movement processes that are constrained by physiological limits and dispersal capability. This can lead to spatial mismatches between local food conditions and prey availability, as in the case of the Baltic Sea, where expansion of hypoxic water has concentrated cod into the remaining suitable habitat that lacks their main prey (Casini et al., 2016). The size of cod in this area has consequently drastically declined, at the same time that cod abundance in the same area has increased (Casini et al., 2016). Thus, spatio-temporal variation in both local habitat conditions and population distributions must be considered when standardizing size-at-age. Given these multiple mechanisms, it can be helpful to attribute variation in annual size-at-age to the various processes contributing to them (e.g. Thorson et al., 2017), in this case, spatial and spatio-temporal variation and local size vs. local abundance.

In this study, we used walleye pollock (Gadus chalcogrammus) in the Bering Sea as a data-rich case study in using a spatio-temporal model to incorporate local-scale spatial and temporal variation in size (i.e. weight) and abundance to estimate population-level size-at-age (i.e. a size-at-age matrix for use in a stock assessment). We evaluate the pattern of size over space and time, whether these patterns are persistent or transient (i.e. varying among years), and their implications for the stock assessment. This fishery is particularly relevant to test a spatio-temporal size-at-age model. Bering Sea walleye pollock is the world’s second-largest single-species fishery by catch weight (FAO, 2021) with an average catch of 1.2 million tons since 1979 and first-wholesale value of $1.55 billion in 2019 (Ianelli et al., 2021). The stock is sampled with extensive fishery-independent surveys, providing a rich data set to evaluate spatio-temporal variability in size-at-age. Finally, size-at-age is an important component of the stock assessment because it converts numbers of fish (which forms the basis of the model dynamics) to biomass and catch recommendations. Although the stock assessment is already transitioning to using spatio-temporal models for estimates of a biomass index and age composition (O'Leary et al., 2020), these models have yet to be applied to size-at-age.

Additionally, climate change is predicted to cause significant shifts in the physical and biological regime (Hunt et al., 2008; Whitehouse et al., 2021) that will greatly impact the spatial distribution of walleye pollock (O'Leary et al., 2020). As the pollock range shifts out of the historical survey footprint, spatial gaps in survey data will increase (O'Leary et al., 2020; O'Leary et al., 2022). In addition to incorporating local-scale variation in growth processes and compensating for local abundance, this spatio-temporal model could therefore also enable better predictions of the pollock stock under future climate change by providing a way to extrapolate to unsampled regions as spatial distributions shift. Further research could explore the impact of including covariates and correlation between ages and hauls, as well as short-term forecasting skill. However, we focus here on the first step (showing how model-based estimators compare with existing methods for expanding size-at-age data) and leave these other topics for future research.

It is often difficult to identify the specific environmental features that cause observed variation in fish populations (Thorson et al., 2017; Dambrine et al., 2021). Because of this, here, we constrained our analysis to quantifying latent spatial and spatio-temporal variation and to retrospective size-at-age data to determine the best analysis protocols, to establish an operating model framework, and to ensure that it could replicate existing estimators and be applied in the stock assessment before delving into covariates and forecasts. Here, we focus on an initial assessment of the pattern of spatial and temporal variation in size-at-age in the Bering Sea pollock, the demographic processes contributing to this variation, and how this model impacts outcomes of the stock assessment.

We seek to answer the following questions:

  1. How does size-at-age vary spatially and temporally, both locally and at the aggregated population-level (i.e. the size-at-age matrix used in the stock assessment)?

  2. How do local spatial and spatio-temporal variation impact the size-at-age matrix?

  3. How does local variation in size vs. abundance contribute to the size-at-age matrix?

  4. How does our spatio-temporal model of the size-at-age matrix affect the stock assessment predictions of population biomass compared to a simple empirical estimate?

To do so, we develop a spatio-temporal model for size-at-age and fish density, and use this to calculate the population-level size-at-age, i.e. the “size-at-age matrix”, that is used as input within many stock-assessment models worldwide (e.g. Kuriyama et al., 2016). We then evaluate the drivers within the model and explore its inclusion in a stock assessment compared to a simple empirical estimate.

Methods

We aimed to estimate annual time series of population size-at-age (i.e. a size-at-age matrix) for the Bering Sea walleye pollock that incorporated local spatial variability in size-at-age, while also compensating for shifts in the relative density of areas with small and large fish. This required an estimate of two response variables—average size and fish density—for each age class for each year in distinct locations across the Bering Sea. Then, to aggregate the local size-at-age at each location up to the population level, we area-expanded fish density to abundance and weighted the size-at-age at each location by the local abundance-at-age in that year to calculate the population-level “size-at-age matrix”. We therefore sought to estimate density and size in a joint model. To accommodate estimating both response variables in a single-joint model, we used a flexible model structure (a Poisson-link delta approach to a generalized linear mixed model) that predicts density-at-age |${n_{a,y,g\;}}$| and size-at-age |${w_{a,y,g}}$| for each year y and age a for distinct locations in space g, while accounting for annual effects and spatial fields that cover spatially and temporally dynamic latent variables. We then evaluated how these estimates might affect key outcomes in a stock assessment model.

Data compilation

The model was fitted to data collected from the NOAA Groundfish Bottom Trawl Survey from 1982 to 2019. Estimating size-at-age relied on the specimen age ai (years), weight vi (grams), and length li (millimeters) measurements that are subsampled as part of hauls in the survey. Surveys were not conducted in 2020, and completed survey data from 2021 were not available at the time of this study. Surveys from 1984 to 1998 did not measure individual fish weights because motion-compensating scales were not widely available at that time (e.g. Wilson and Armestead, 1991; Walters, 1997; Nebenzahl and Goddard, 2000). After 1991, specimen weights were sampled sporadically, and were consistently measured from 1999 onwards. For 1984–1999, weight was estimated from measured lengths using sex-specific weight–length conversion parameters (Ianelli et al., 2021). Individual fish from the northern Bering Sea survey were not aged by the time this study was completed. Additionally, there were no specimen data available for age-14 in 2018 and age-15+ in 1986 and 2019. The density-at-age data that we used are the age-stratified, density-dependence-corrected catch-per-unit-effort (CPUE) per haul (numbers of fish hectare−1) ci [see Kotwicki et al. (2014) for details]. Figure 1 shows the spatial distribution of survey data for a subset of years, and the general geographic extent of CPUE and specimen weight or length data are similar. See Supplementary Figure 1 for full time series of locations of specimen weight and length data and Supplementary Figure 2 for CPUE data.

Bottom trawl survey locations in the Bering Sea for subset of years to show extent of survey over time.
Figure 1.

Bottom trawl survey locations in the Bering Sea for subset of years to show extent of survey over time.

Model structure

Average local size |${w_{g,a,y}}$| for each year y and age a for distinct locations in space g was estimated using a log-linked generalized linear mixed model. The model includes a single linear predictor that contains terms for temporal variation |$\beta $|⁠, spatial variation |$\omega $|⁠, and spatio-temporal variation |$\varepsilon $|⁠.
(1)
Temporal variation was estimated as a fixed effect for each age a and year y, and represents a fixed intercept for each age and year. Spatial and spatio-temporal variation were estimated as Gaussian random effects with mean 0 and were assumed to follow a multivariate normal distribution with covariance matrices |$\sigma _{\rm{\omega }}^2{{\bf{R}}_{\bf{\omega }}}$| and |$\sigma _\varepsilon ^2{{\bf{R}}_\varepsilon }$|⁠. Covariance matrices include the correlation between locations, calculated using a Matérn function in the correlation matrices |${{\bf{R}}_{\bf{\omega }}}$| and |${{\bf{R}}_\varepsilon }$|⁠, which are governed by a shared decorrelation rate |${\rm{\kappa }}$| (i.e. the distance that locations are uncorrelated), and the estimated marginal variances |${\sigma _\omega }$| and |${\sigma _\varepsilon }.\;$| Spatio-temporal variation also follows an auto-regressive process, governed by the autocorrelation parameter |${\rho _\varepsilon }$|⁠.
(2)
(3)
The measured size |${v_i}$| is assumed to follow a lognormal distribution, with log-mean |$\log ( {{w_{{g_i},{a_i},{y_i}}}} )$| and log-variance |$\sigma _{\rm{v}}^2$|⁠.
(4)
In order to weight the local size estimates from this size-at-age model by the local abundance at each location in each year, we additionally needed a similar spatio-temporal model of density-at-age |${n_{g,a,y\;}}$| for each year y and age a for distinct locations in space g. However, catch rates contain many zeroes for certain ages in certain years (Martin et al., 2005). To accommodate for zero-inflation in the model of abundance, we used a Poisson-link delta approach to a generalized linear mixed model (Thorson, 2018). This approach calculates the probability of encounter |${r_{1,i}}$| from a Poisson distribution, whose rate parameters are governed by linear predictor |${p_{1,i}}$| (i.e. a complementary log–log link function), and the density if encountered |${r_{2,i}}$| from the Poisson distribution modeled as |${p_{1,i}}$|⁠, and a gamma-distributed multiplier of counts, whose mean and variance are determined by a linear predictor |${p_{2,i}}.$| Thus, the encounter rate |${r_{1,i}}$| and density if encountered |${r_{2,i}}$| equal:
(5)
(6)
Final estimated density |${n_{i\;}}$| is then the product of |${r_{1,i}}$| and |${r_{2,i}}$|⁠, such that |$\log ( {{n_i}} ) = {p_{1,i}} + {p_{2,i}}$|⁠.
Because both linear predictors are included in the calculation of positive catch rates |${r_{2,i}}$|⁠, including the three variation terms in the first linear predictor includes them in both |${r_{1,i}}$| and |${r_{2,i}}$|⁠. To simplify the model for computational efficiency, temporal, spatial, and spatio-temporal variation are included in the first linear predictor, |${p_{1,i}}$| , while only temporal variation was included in the second linear predictor. When using the Poisson-link delta model, the variance associated with the second linear predictor is typically small, and we can assume that excluding the terms in the second linear predictor will have negligible impacts on the model estimates (see Thorson 2018).
(7)
(8)
The probability of not encountering an age-class in a given sample [i.e. Pr(⁠|${C_i}$|=0)] equals |$1 - {r_{1,i}}$|⁠, and the positive catch rates |${C_i}$| > 0 were assumed to follow a gamma distribution, governed by shape parameter |$k$| and scale parameter |$\theta $|⁠, with mean |${r_{2,i}}$| and coefficient of variation (CV) |${\sigma _r}$|⁠, where |$k = \frac{1}{{\sigma _r^2}}$| and |$\theta $|=|${r_{2,i}}\; \times \;\sigma _r^2$|⁠.
(9)

We implemented this joint spatio-temporal size- and density-at-age model using the R package VAST version 3.8.2 (Thorson, 2019). This package allows the flexibility to do a multivariate spatio-temporal model, essentially treating each size- and density-at-age as a different response variable (i.e. category). Each category is specified with its probability distributions and link function, and separate spatial (i.e. omega term) and spatio-temporal (i.e. epsilon term) variation are estimated for each size- and density-at-age. Spatial and spatio-temporal terms are not correlated across ages. This package also allowed the flexibility to estimate both size- and density-at-age by using a Poisson-link delta approach that estimated density-at-age as in Equations (5) to (9), but essentially fixed the encounter probability for size-at-ages at 1 (i.e. 100% encountered) to ultimately estimate size-at-age by a single linear predictor. For additional details on how this joint model was coded in VAST, see Supplementary Methods Section 2. Additional information on code and documentation of the VAST package is available at https://github.com/James-Thorson-NOAA/VAST.

Ultimately, the model extrapolated size and abundance by predicting density-at-age |${n_{g,a,y\;}}$| and size-at-age |${w_{g,a,y}}$| for each year y and age a for distinct locations in space g. The model estimated spatial and spatio-temporal effects across 500 points (i.e. termed “knots” in VAST) that were uniformly distributed across the Bering Sea and used bilinear interpolation to cover the entire geographic extent of the model (in this case, the Bering Sea) between these points with 51769 distinct locations (i.e. grid cells g). We established that 500 knots were sufficient by qualitatively comparing model predictions over a range of knots. Density- and size-at-age were estimated for each of the years 1982–2019 at each of these 51769 locations spanning the region. Approximating the multivariate normal distribution was accomplished using the SPDE method (Lindgren et al., 2011), and the matrices involved were calculated using R-INLA (Lindgren and Rue, 2015). See Supplementary Figure 3 for the map of the SPDE mesh, Supplementary Table 2 for estimated model parameters, and Table 1 for a summary of definitions of variable names and symbols.

Table 1.

Variable names and symbols used in data processing and model outputs.

NameSymbolUnits
Indices
Specimen (data)i
Haul (data)h
Age (data and model)aYear
Year (data and model)y
Grid cell (model)g
Data inputs
Fish specimen size (weight)viGrams
Fish specimen fork lengthliMillimeters
Catch rate (density-corrected CPUE)ciNumber of fish hectare−1
Model components
Temporal variation|$\beta $|
Spatial variation|$\omega $|
Spatio-temporal variation|$\varepsilon $|
Linear predictor|$p$|
Encounter probability|${r_1}$|
Density if encountered|${r_2}$|
Autocorrelation for spatio-temporal variation|${\rho _\varepsilon }$|
Spatial correlation matrix for spatio-temporal variation|${{\bf{R}}_\varepsilon }$|
Spatial correlation matrix for spatio-temporal variation|${{\bf{R}}_{\bf{\omega }}}$|
Shape parameter of gamma distribution|$k$|
Scale parameter of gamma distribution|$\theta $|
Model outputs
Fish density (i.e. density)nNumber of fish hectare−1 (or km−2)
Abundance|$\hat n$|Number of fish
Size (i.e. weight)wGrams
Abundance-weighted size|$\hat w$|Kilograms
AreaAKilometer2
NameSymbolUnits
Indices
Specimen (data)i
Haul (data)h
Age (data and model)aYear
Year (data and model)y
Grid cell (model)g
Data inputs
Fish specimen size (weight)viGrams
Fish specimen fork lengthliMillimeters
Catch rate (density-corrected CPUE)ciNumber of fish hectare−1
Model components
Temporal variation|$\beta $|
Spatial variation|$\omega $|
Spatio-temporal variation|$\varepsilon $|
Linear predictor|$p$|
Encounter probability|${r_1}$|
Density if encountered|${r_2}$|
Autocorrelation for spatio-temporal variation|${\rho _\varepsilon }$|
Spatial correlation matrix for spatio-temporal variation|${{\bf{R}}_\varepsilon }$|
Spatial correlation matrix for spatio-temporal variation|${{\bf{R}}_{\bf{\omega }}}$|
Shape parameter of gamma distribution|$k$|
Scale parameter of gamma distribution|$\theta $|
Model outputs
Fish density (i.e. density)nNumber of fish hectare−1 (or km−2)
Abundance|$\hat n$|Number of fish
Size (i.e. weight)wGrams
Abundance-weighted size|$\hat w$|Kilograms
AreaAKilometer2
Table 1.

Variable names and symbols used in data processing and model outputs.

NameSymbolUnits
Indices
Specimen (data)i
Haul (data)h
Age (data and model)aYear
Year (data and model)y
Grid cell (model)g
Data inputs
Fish specimen size (weight)viGrams
Fish specimen fork lengthliMillimeters
Catch rate (density-corrected CPUE)ciNumber of fish hectare−1
Model components
Temporal variation|$\beta $|
Spatial variation|$\omega $|
Spatio-temporal variation|$\varepsilon $|
Linear predictor|$p$|
Encounter probability|${r_1}$|
Density if encountered|${r_2}$|
Autocorrelation for spatio-temporal variation|${\rho _\varepsilon }$|
Spatial correlation matrix for spatio-temporal variation|${{\bf{R}}_\varepsilon }$|
Spatial correlation matrix for spatio-temporal variation|${{\bf{R}}_{\bf{\omega }}}$|
Shape parameter of gamma distribution|$k$|
Scale parameter of gamma distribution|$\theta $|
Model outputs
Fish density (i.e. density)nNumber of fish hectare−1 (or km−2)
Abundance|$\hat n$|Number of fish
Size (i.e. weight)wGrams
Abundance-weighted size|$\hat w$|Kilograms
AreaAKilometer2
NameSymbolUnits
Indices
Specimen (data)i
Haul (data)h
Age (data and model)aYear
Year (data and model)y
Grid cell (model)g
Data inputs
Fish specimen size (weight)viGrams
Fish specimen fork lengthliMillimeters
Catch rate (density-corrected CPUE)ciNumber of fish hectare−1
Model components
Temporal variation|$\beta $|
Spatial variation|$\omega $|
Spatio-temporal variation|$\varepsilon $|
Linear predictor|$p$|
Encounter probability|${r_1}$|
Density if encountered|${r_2}$|
Autocorrelation for spatio-temporal variation|${\rho _\varepsilon }$|
Spatial correlation matrix for spatio-temporal variation|${{\bf{R}}_\varepsilon }$|
Spatial correlation matrix for spatio-temporal variation|${{\bf{R}}_{\bf{\omega }}}$|
Shape parameter of gamma distribution|$k$|
Scale parameter of gamma distribution|$\theta $|
Model outputs
Fish density (i.e. density)nNumber of fish hectare−1 (or km−2)
Abundance|$\hat n$|Number of fish
Size (i.e. weight)wGrams
Abundance-weighted size|$\hat w$|Kilograms
AreaAKilometer2

Annual population size-at-age

From these 51769 local estimates of size-at-age for each year, we wanted an average annual population-level size-at-age (i.e. “size-at-age matrix”) that incorporated local variation in size and compensated for the local abundance. This is because the stock assessment model requires a size-at-age matrix with dimensions age × year. To that end, for each year, we calculated a weighted average of size-at-age, weighted by the local abundance. First, the modelled density |${n_{g,a,y\;}}$| per age a, year y, and grid cell g (in number of fish hectare−1) was expanded by the area |${A_g}$| of the grid cell (in kilometers2) to calculate area-expanded abundance |${\hat n_{g,a,y}}$| (with units numbers of fish) using Equation (10). |${\hat n_{g,a,y}}$| was then summed across all grids per year and per age to calculate the annual, population-level area-expanded abundance |${\hat n_{a,y}}$|⁠.
(10)
Second, a weighted average of size for each age class and year, |${\hat w_{a,y}},$| was calculated by weighting size-at-age |${w_{a,y,g}}$| by the abundance associated with this extrapolation point |${n_{g,a,y\;}} \times {A_g}$|⁠.
(11)

Question 1: spatial and temporal patterns of size-at-age

We descriptively evaluated spatial and temporal patterns in the model predictions of size and abundance at the local level and in population-level size-at-age matrix. We also calculated the CV in size for each age class to quantify the extent of variation and compare patterns in size-at-age.

We were also interested in comparing our spatio-temporal model of the size-at-age matrix to a non-spatially explicit simple estimate. For this simple empirical estimate, we calculated the weighted average of size-at-age per year from the observed size data, weighting mean size-at-age in each haul by catch-at-age in that haul from the survey data. (See Supplementary Material for detailed equations.)

Question 2: impact of spatial and spatio-temporal variation on size-at-age matrix

We sought to compare how much spatial vs. spatio-temporal variation impacted the model predictions of the size-at-age matrix used in the Bering Sea walleye pollock stock assessment. A greater contribution of spatio-temporal variation indicates that there are important transient annual changes in habitat suitability. To do so, we used a method similar to adjusted predictions or marginal effects (e.g. Williams, 2012; Bornmann and Williams, 2013; Mize, 2019) and to the counterfactual approach in Thorson et al. (2017), where model predictions are compared when some variables in the model are held constant, while the variable of interest fluctuates as a way to evaluate the effect of each variable within the model. In our study, the spatial term |${\omega _{g,a\;\;}}$| and spatio-temporal term |${\varepsilon _{g,y,a\;\;}}$| from each grid cell was each alternately replaced with the respective mean to calculate a reduced model (i.e. no spatial or no spatio-temporal variation term) of population size-at-age |${\hat w_{a,y}}$| for the size-at-age matrix. Because spatial and spatio-temporal terms are random effects with means of 0, this essentially “zeroes out” the term and allows us to see from the change in model prediction the relative impact of spatial or spatio-temporal variation in the model. We used this approach, rather than compare predictions between alternative model fits, because we aimed to compare the relative strength of each type of variation within our model structure when both spatial and spatio-temporal variation are estimated. Our approach allows us to quantify the effect of each model component when used in the inverse-link transformed linear predictor, and hence, provides a more useful summary that simply interpreting the estimated standard deviation of each model component as its “importance”. The % difference between |${\hat w_{a,y}}$| in the reduced model (i.e. no-spatial or no-spatio-temporal variation) from the full model (i.e. all forms of variation) was then calculated for each year. We interpreted the greater % difference between the simplified and full model as indicating a greater impact of that source of variation.

Question 3: contribution of local variation in size vs. abundance to size-at-age matrix

We similarly evaluated how the size-at-age matrix was impacted by accounting for the local variation in size and by local abundance. To do so, we alternately replaced the grid-level abundance-at-age and grid-level size-at-age with their average across all years, and re-calculated the size-at-age matrix using these values. Specifically:

  • To isolate the impact of correcting for local abundance, the local area-expanded abundance |${\hat n_{g,a,y}}$| was replaced with the mean abundance-at-age for each grid cell across all years, |${\overline {\hat n} _{g,a}}$|⁠.

  • Alternately, to isolate the effect of local size-at-age, the local size, |${w_{g,a,y}}$|⁠, was replaced with the average size for each grid cell for each age class |${\bar w_{{\rm{g}},{\rm{\;a}}}}$| across years and subsequently used to calculate annual size-at-age |${\hat w_{a,y}}$|⁠.

The % difference between |${\hat w_{a,y}}$| in the reduced model (i.e. no local temporal variation in size or no local temporal variation in density) from the full model was calculated for each year.

Question 4: stock assessment comparison

We evaluated how using our spatio-temporal model to calculate the size-at-age matrix |${\hat w_{a,y}}$|—incorporating local spatio-temporal variation in size and weighting by abundance—impacts outcomes of the stock assessment. We did so by comparing the walleye pollock stock assessment model (see Ianelli et al., 2021) using our spatio-temporal size-at-age matrix to using the simple empirical estimate (see above). Because there was no specimen data available for age-14 in 2018 and age-15+ in 1986 and 2019, we used average values since there were no model estimates for those age classes in those years. We ran the Bering Sea walleye pollock stock assessment model using each size-at-age matrix method and compared the total estimated biomass and the predicted numbers of fish in the stock from the stock assessment model. We compared the fit of the stock assessment model using each size-at-age method to the trawl survey biomass using Akaike information criteria (AIC), where a value of <2 is weak, >10 is very strong, and between these is intermediate (Burnham and Anderson, 2002) .

Results

Question 1: spatial and temporal patterns of size-at-age

Size-at-age varied across the Bering Sea and among years. Across the Bering Sea, we found some persistent spatial patterns (i.e. spatial variation). For instance, age-1 pollock size reached 30–40 g in the western extent of its range, while averaging only 10–20 g in the eastern extent of the range (Figure 2). Age-9 pollock were typically larger at the eastern end of the range than the western, reaching 2000 g in some years in the eastern end, while less than half that size in the western end (Figure 2). However, the spatial patterns of size differed among years (i.e. spatio-temporal variation). For example, an area of particularly higher age-1 size appeared in the northeast in 2019 (Figure 2), and some years exhibited more variability in size across the region. For example, while there was great spatial variation and stratification of size in age-9 pollock in 1982 and 2000, there was a more uniform size-at-age across the entire Bering Sea in 1990 and 2019.

Model estimates of area-expanded abundance log(abundance, numbers of fish) and weight (grams) for age-1 and age-9 walleye pollock (G. chalcogrammus) for subset of years (1982, 1990, 2000, 2010, and 2019).
Figure 2.

Model estimates of area-expanded abundance log(abundance, numbers of fish) and weight (grams) for age-1 and age-9 walleye pollock (G. chalcogrammus) for subset of years (1982, 1990, 2000, 2010, and 2019).

The overall magnitude of spatial variation in size also differed among age classes. The youngest age classes showed the greatest variation across grid cells (Figure 3), with a CV of 42% for age-1 and 39% for age-2 pollock. There was a marked decline in variation for older ages. CV dropped to 30% for age-3 and further declined to ∼18–22% for each age classes 4–15+ .

Spatial CV of estimated size-at-age of walleye pollock (G. chalcogrammus) in the Bering Sea for each age class, calculated as the standard deviation across all grid cells and all years divided by the mean across grid cells and years.
Figure 3.

Spatial CV of estimated size-at-age of walleye pollock (G. chalcogrammus) in the Bering Sea for each age class, calculated as the standard deviation across all grid cells and all years divided by the mean across grid cells and years.

There was also a trend of declining size-at-age over the study period for all but the youngest age classes at both the local and population level (i.e. the size-at-age matrix). Beginning in 2015 through 2019, there was a distinct, uniform decline in weight for age-8 pollock and older across the entire Bering Sea, while this was not evident in younger ages (e.g. see Supplementary Figure 4). At the population level, size-at-age from the first modelled year (1982) to the last modelled year (2019) declined by 24% for age-1 pollock, increased by 12–60% for age-2 through age-5 pollock, and declined by 15–30% for age classes 6–8, and 35–50% for age classes 9 and older (e.g. Figure 4 and Supplementary Figure 5).

Size-at-age matrix (i.e. annual abundance-expanded population size-at-age) from model estimates (dashed line) and standard error (grey-shaded), compared to the simple empirical estimate (i.e. the simple empirical average, a non-spatially explicit weighted average of observations) (solid black line) for ages-1 (top panel) and 9 (bottom panel) walleye pollock. See Supplementary Figure 4 for all age classes.
Figure 4.

Size-at-age matrix (i.e. annual abundance-expanded population size-at-age) from model estimates (dashed line) and standard error (grey-shaded), compared to the simple empirical estimate (i.e. the simple empirical average, a non-spatially explicit weighted average of observations) (solid black line) for ages-1 (top panel) and 9 (bottom panel) walleye pollock. See Supplementary Figure 4 for all age classes.

Our spatio-temporal size-at-age matrix differed from the simple empirical estimate (i.e. the non-spatially explicit, simple empirical average) of size-at-age. Though this non-spatial estimate somewhat tracked fluctuations in the model size-at-age, it often fell outside the confidence range (Figure 4 and Supplementary Figure 5).

Question 2: impact of spatial and spatio-temporal variation on size-at-age matrix

The size-at-age matrix was impacted by both spatial and spatio-temporal variability. Both latent persistent spatial and transient spatio-temporal patterns in size across the Bering Sea contributed to the fluctuations in the population-level size-at-age matrix. Population-level size-at-age that excluded spatial or spatio-temporal variation differed from the full model (Figure 5 and Supplementary Figure 6), though the magnitude of the impact varied between age classes. Spatial variation had the largest impacts in age-1 and ages 5–7, with 10–25% average impact on annual size-at-age. However, for the remaining age classes, the impact of spatial variation was relatively minor, with an average difference of 8–14% for age-8–14, and almost zero impact on age-3 and age-15. Spatio-temporal variation had an overall lower impact on more age classes, ranging from 2 to 8% for age-1 and age-4–15, and only had a larger impact on age-2 and age-3 pollock (14 and 20%, respectively). Overall, spatial and spatio-temporal variation had relatively similar magnitude of impact on the size-at-age matrix.

The size-at-age matrix (i.e. annual abundance-expanded population size-at-age) from the full model, including spatial, spatio-temporal, and temporal variation (solid black line) compared to (1) removing spatial variation (dashed line) and (2) removing spatio-temporal variation (dotted line). See Supplementary Figure 5 for all age classes.
Figure 5.

The size-at-age matrix (i.e. annual abundance-expanded population size-at-age) from the full model, including spatial, spatio-temporal, and temporal variation (solid black line) compared to (1) removing spatial variation (dashed line) and (2) removing spatio-temporal variation (dotted line). See Supplementary Figure 5 for all age classes.

Question 3: contribution of local variation in size vs. abundance to size-at-age matrix

Annual fluctuations in the weight-at-age matrix were almost entirely due to changes in local size between years, rather than shifts in local fish abundance (Figure 6). Inter-annual fluctuations in size-at-age were markedly dampened when local variation in size-at-age was removed, but were not substantially changed when local variation in abundance was removed. Incorporating local variation in abundance-at-age but holding local size constant caused an 11–19% average difference (Figure 6 and Supplementary Figure 7). Including local changes in fish size while holding local abundance-at-age constant on average caused only a 2–8% average difference (Figure 6). Overall, accounting for the changes in local size-at-age processes year-to-year had a substantial impact when extrapolating to the population level, while compensating for local changes in fish density did not greatly change the size-at-age matrix value.

The population size-at-age matrix from the full model (i.e. calculating the abundance-expanded size-at-age with both local size-at-age and local abundance-at-age) (solid black line) compared to the size-at-age matrix calculated using (1) local abundance but mean size-at-age (i.e. showing the contribution of local size-at-age) (dotted line) and (2) mean abundance but local size-at-age (i.e. showing the contribution of local abundance) (dashed line). See Supplementary Figure 6 for all age classes.
Figure 6.

The population size-at-age matrix from the full model (i.e. calculating the abundance-expanded size-at-age with both local size-at-age and local abundance-at-age) (solid black line) compared to the size-at-age matrix calculated using (1) local abundance but mean size-at-age (i.e. showing the contribution of local size-at-age) (dotted line) and (2) mean abundance but local size-at-age (i.e. showing the contribution of local abundance) (dashed line). See Supplementary Figure 6 for all age classes.

Question 4: stock assessment comparison

The spatio-temporal size-at-age matrix did marginally change estimated annual biomass and numbers of fish from the walleye pollock stock assessment (Figure 7) compared to the simple empirical mean. The stock assessment model using our spatio-temporal size-at-age resulted in a slightly improved fit than a model using the simple empirical estimate (ΔAIC = 3.6). In particular, the spatio-temporal model resulted in marginally better predicted biomass in the year 2003 (Figure 7), when the observed biomass in the bottom trawl survey was anomalously higher than the model estimates for that year. The spatio-temporal model was able better avoid the underestimation and capture the higher biomass that was a closer fit to the survey observation. A similar pattern is seen in 1991, though there is less divergence between the model and observation than 2003. The spatio-temporal size-at-age matrix generally estimated lower size-at-age than the simple empirical estimate in 1991, except for much higher estimates in the oldest age-classes of 13–15+ . In 2003, size-at-age was also higher in the spatially explicit model in the oldest age classes, but also in some younger ages (Supplementary Figure 5). There was essentially no difference in the estimated numbers of fish in the stock by the assessment model between size-at-age methods until the most recent several years (2015–2021), when the numbers of fish estimated using the spatio-temporal size-at-age are consistently higher than the simple empirical estimate (Figure 4). This suggests that local-scale variation in weight did manifest in the demographic patterns of the population, and that accounting for these structures scaled up to impact the estimation of population-level metrics.

Walleye pollock biomass from bottom trawl survey (top) and predicted numbers of fish (bottom) estimated from stock assessment model (Ianelli et al., 2021) using spatio-temporal size-at-age matrix (blue) and from the simple empirical estimate (a non-spatially explicit weighted average of observations) (red), compared to observed biomasses (black points).
Figure 7.

Walleye pollock biomass from bottom trawl survey (top) and predicted numbers of fish (bottom) estimated from stock assessment model (Ianelli et al., 2021) using spatio-temporal size-at-age matrix (blue) and from the simple empirical estimate (a non-spatially explicit weighted average of observations) (red), compared to observed biomasses (black points).

Discussion

We developed and demonstrated a population-level size-at-age matrix that accounts for both inherent differences in size-at-age between habitats across the region and movement of fish (i.e. shifts in distribution) among patches with smaller or larger size. Local spatial variation in size-at-age, rather than local changes in abundance, largely governed variability in the resulting weight-at-age matrix. Accounting for spatial and spatio-temporal variation in local size-at-age therefore considerably impacts the estimated population size-at-age. Spatio-temporal models, such as with VAST as done here, of abundance (e.g. Fenske et al., 2020) and age composition (e.g. Thorson and Haltuch, 2018; O'Leary et al., 2020) have been used in stock assessments previously, but to our knowledge, this is the first example of the approach being applied to size-at-age. While the stock assessment in our study was relatively robust to the method used to calculate survey size-at-age, climate change will make it increasingly important for fisheries to account for the rapid pace of distribution shifts and environmental changes that species are likely to experience. The overall framework implemented here (spatio-temporal estimation and abundance-expansion of size-at-age) could be widely applied to other systems and stock assessments to evaluate the sensitivity of estimates to changing species distributions and local changes in productivity and growth.

This study highlights large decreases in size-at-age for pollock. Our model showed a general trend of a decline in size-at-age between the beginning and end of the study period in all but the youngest age classes of pollock. While age-2 through age-5 pollock did not see a decline, sizes of age-1 and ages-6–15+ in the final year of the study period (2019) were 15–50% lower than at the beginning (1982). This aligns with recent reports from fishers, who noted smaller than usual pollock in 2020 and 2021 catches (Ianelli et al., 2021).

This decline may be due to changes in environmental conditions or size-selective fishing. Declines in fish growth and size with ocean warming are expected due to increased metabolic demand and reduced oxygen availability in higher sea temperatures (Daufresene et al., 2009; Gardner et al., 2011) and have been widely documented (e.g. Baudron et al., 2014; van Rijn et al., 2017). For instance, Pacific halibut mean size-at-age has declined at a similar magnitude to our study (Holsman et al., 2018), and has been attributed to both sea temperature (Holsman et al., 2018) as well as interspecific competition, density-dependent effects, and size-selective fishing (Clark and Hare, 2002; Sullivan et al., 2018). Higher temperatures are also likely to reduce body size in larger and older individuals more than smaller and younger individuals (Ikpewe et al., 2020; Lindmark et al., 2022; Oke et al., 2022). Over the last several years (∼2016–2021), the Bering Sea has experienced warmer than average temperatures (Ianelli et al., 2021). The age-class patterns we observed and recent environmental conditions of the Bering Sea would be consistent with an age-specific growth response to environmental conditions.

Size-selective fishing may also have contributed to the decline in size-at-age. Fishery-induced declines in size-at-age have been documented in numerous fisheries (Sharpe and Hendry, 2009) such as Pacific halibut (Sullivan, 2016, Sullivan et al., 2018) and northern cod (Krohn and Kerr, 1997). Size-selective fishing can lead to smaller body size through evolutionary pressure (e.g. Swain et al., 2007; Jorgensen et al., 2009) causing heritable genetic change (Uusi-Heikkila et al., 2015), phenotypic plasticity (Fenberg and Roy, 2008), or the higher removal of faster-growing individuals of a population (i.e. the “Rosa Lee phenomenon”, Kraak et al., 2019). In the Bering Sea pollock fishery, vessel operators pursue pollock that are optimal for cost-efficient processing and market demand, and some vessels have modified gear to exclude smaller pollock (Ianelli et al., 2021). Further exploration of the links between spatio-temporal variation in size-at-age of walleye pollock in the Bering Sea and environmental conditions would help clarify the possible cause of declining size-at-age.

This spatio-temporal model of size-at-age can improve stock assessments in several ways. Our estimated weight-at-age matrix incorporated the extensive local variation in size across the eastern Bering Sea. By jointly estimating size and abundance at a fine spatio-temporal scale, we also compensated for uneven distributions in fish abundance and shifts in distribution of fish throughout the region when expanding the local size estimates to the population size-at-age matrix. While the stock assessment was not substantially impacted by the method for generating the size-at-age matrix, the spatio-temporal and abundance-expanded size-at-age matrix did slightly improve the fit of the model to the biomass survey data. For instance, the spatio-temporal size-at-age was able to slightly improve the model’s ability to capture the anomalously high biomasses observed in the bottom trawl survey in 2003 (and to a lesser extent in 1991). Lastly, a spatio-temporal model, as done in this study, can help extrapolate in time and space to account for shifting survey and stock extents (O'Leary et al., 2020, 2022; Adams et al., 2021). For instance, in this study, weight observations of pollock in the northern Bering Sea were not available, and there is sparse data from that region because of limited surveys. Additionally, our joint modelling approach allows the size-at-age matrix to include both uncertainty in estimates of size and uncertainty in estimates of abundance. Spatio-temporal approaches have become more common practice in some components of the stock assessment (such as biomass and age composition). While we did not see a marked change in the stock assessment outcome from our spatio-temporal method to size-at-age, accounting for spatio-temporal variation in other demographic metrics, such as size-at-age in this study, are the logical next step and will improve methodological consistency across a stock assessment.

Extending our spatio-temporal framework of size-at-age to fishery size data could be particularly beneficial because fishery-dependent data tend to have more pronounced spatial variation in fishing effort. There are differences in selectivity between fleets (e.g. gear characteristics such as depth of fishing) and unbalanced sampling in areas of higher, or lower, fish abundance due to fishers preferentially targeting areas based on convenience, safety, profitability, etc. (Maunder et al., 2020). A spatio-temporal model of size-at-age from fishery data consequently would require additionally accounting for catchability and selectivity when weighting size-at-age by catch (Maunder et al., 2020; Thorson et al., 2020). Stock assessments can be sensitive to how fishery size-at-age is estimated (Punt and Smith, 2001). An accurate fishery size-at-age is necessary for correctly determining the number of fish caught by the fishery when landings are measured in weight and, consequently, how actual fishing mortality compares to management limits (Ianelli et al., 2021). This study shows that it is feasible to take the grid-level output of size and abundance estimates from VAST and calculate a population average size-at-age, first weighting density by area and then weighting size by abundance. Including catchability covariates and weighting size by catch and fitting to fishery data would be a tractable extension of the model used in this study.

Our model showed extensive local-scale differences in size-at-age across the region and across years, suggesting that there is significant variation in local growth processes. Certain areas may have intrinsic conditions that make them more or less favourable for growth (e.g. Williams et al., 2003, Begg and Martensdottir, 2002), contributing to the persistent spatial differences seen. Spatial patterns in pollock size have been previously observed; for instance, pollock in the northwest area have typically been found to be smaller than those in the southeast (Ianelli et al., 2021). Additionally, we saw annual differences in the spatial patterns that may indicate more transient changes in habitat conditions. We also saw some evidence of local shifts in abundance, though this was not as predominant as local variation in size-at-age. Passive dispersal and habitat selection can impact size-at-age if fish move to optimize foraging resources, to access spawning grounds, or to avoid competition or predators (Hanselman et al., 2015). Fish may also make dynamic migratory decisions based on energetic status, resulting in larger individuals migrating to different habitats than smaller individuals, such as has been seen in sablefish (Hanselman et al., 2015). In this study, however, local changes in growth processes, rather than shifts in distribution, seemed to be the principal driver of size-at-age.

Our spatio-temporal model of size-at-age also informs our understanding of how these local spatial and temporal patterns in size affect population-scale productivity. This has been identified as a “Grand Habitat Challenge” (Thorson et al.,2021). Here, we developed a model of local size at a fine spatial scale, and expanded it to population-level abundance and productivity through the stock assessment. We show that persistent habitat conditions (spatial variation) and annual fluctuations in the suitability of a habitat for growth (spatio-temporal variation) caused local-scale variation in size-at-age over space and between years. Incorporating the variation in local size-at-age was the main driver of the annual population size-at-age (i.e. the size-at-age matrix). This suggests that variation in local growth processes across years did have an impact on population-level demographic metrics. In this study, we therefore show a feasible approach to incorporate fine-scale local demographic processes into population-level assessments.

Being able to predict how demographic rates, such as growth and size, may respond to environmental change will also be crucial for incorporating the effects of climate change into stock assessment models (Punt et al., 2021). In the Bering Sea, climate change is likely to alter environmental conditions known to impact walleye pollock demographics. Declining sea ice cover and associated declines in zooplankton may reduce food availability and impact growth processes (Hunt et al., 2008; Mueter and Litzow, 2008; Hunt et al., 2022) and recruitment (Mueter et al., 2011), and reduced bottom sea temperatures (i.e. the “cold pool extent) impacts a suite of demographic characteristics (Gruss et al., 2021). A model that can estimate size-at-age at fine spatial and temporal scales will enable better predictions of size-at-age, and the associated impacts on yield-per-recruit and other elements of the stock assessment, under future climate scenarios. For simplicity, in this study, we did not integrate environmental or catchability covariates or forecast future size-at-age. Rather, we focused on an initial assessment of the feasibility of applying spatio-temporal models to size-at-age and abundance for stock assessments. It can often be difficult to identify the specific environmental features that cause observed variation in fish populations (Thorson et al., 2017; Dambrine et al., 2021), and addressing the latent spatio-temporal as done here was a tangible first step. The framework presented here can be expanded to incorporate environmental covariates and used for short-term forecasting.

Supplementary data

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

ACKNOWLEDGEMENTS

We would like to specifically thank Cecilia O'Leary and Caitlin Allen-Akelsrud for their help and guidance on pulling together the various data sources. Special thanks to the bottom trawl survey team and Age and Growth Program at the NOAA Alaska Fisheries Science Center for the data that make this work possible. We also thank Paul Spencer and Christine Stawitz and the two anonymous reviewers for helpful comments on a previous draft.

Funding

This publication was partially funded by the Cooperative Institute for Climate, Ocean, & Ecosystem Studies (CIOCES) under NOAA Cooperative Agreement NA20OAR4320271, Contribution No. 2022-1205 and by the National Pacific Research Board Grant #1805. J.I. was supported partly by a CICOES Graduate Student Award (Contribution No. 2022-1205) and by the National Pacific Research Board Grant #1805.

Data availability

The data underlying this article will be shared on reasonable request to the corresponding author.

Conflict of interest

The authors have no conflicts of interest to declare.

References

Adams
C.
,
Brooks
E.
,
Legault
C.
,
Barrett
M.
,
Chevrier
D.
2021
.
Quota allocation for stocks that span multiple management zones: analysis with a vector autoregressive spatiotemporal model
.
Fisheries Management and Ecology
,
28
:
417
427
.

Baudron
A.
,
Needle
C.
,
Rijnsdorp
A.
,
Marshall
C.
2014
.
Warming temperatures and smaller body sizes: synchronous changes in growth of North Sea fishes
.
Global Change Biology
,
20
:
1023
1031
.

Begg
G.
,
Marteinsdottir
G.
2002
.
Environmental and stock effects on spatial distribution and abundance of mature cod Gadus morhua
.
Marine Ecology Progress Series
,
229
:
345
262
.

Bornmann
L.
,
Williams
R.
2013
.
How to calculate the practical significance of citation impact differences? An empirical example from evaluative institutional bibliometrics using adjusted predictions and marginal effects
.
Journal of Informetrics
,
7
:
562
574
.

Breivik
O.
,
Aanes
F.
,
Sovik
G.
,
Aglen
A.
,
Mehl
S.
,
Johnsen
E.
2021
.
Predicting abundance indices in areas without coverage with a latent spatio-temporal Gaussian model
.
ICES Journal of Marine Science
,
78
:
2031
2042
.

Burnham
K.
,
Anderson
D.
2002
.
Model Selection and Inference
.
Springer-Verlag
,
New York, NY
.

Casini
M.
,
Kall
F.
,
Hansson
M.
,
Plikshs
M.
,
Baranova
T.
,
Karlsson
O.
,
Lundstrom
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
.

Clark
W.
,
Hare
S.
2002
.
Effects of climate and stock size on recruitment and growth of Pacific halibut
.
North American Journal of Fisheries Management
,
22
:
852
862
.

Dambrine
C.
,
Woillez
M.
,
Huret
M.
,
de Pontual
H.
2021
.
Characterising essential fish habitat using spatio-temporal analysis of fishery data: a case study of the European seabass spawning areas
.
Fisheries Oceanography
,
30
:
413
428
.

Daufresne
M.
,
Lengfellner
K.
,
Sommer
U.
2009
.
Global warming benefits the small in aquatic ecosystems
.
Proceedings of the National Academy of Sciences
,
106
:
12788
12793
.

DeVries
D.R.
,
Frie
R.V.
1996
.
Determination of age and growth
. In
Fisheries Techniques
, 2nd edn, pp.
483
512
.. Ed. by
B.R.
Murphy
,
D.W.
Willis
.
American Fisheries Society
,
Bethesda, MD
.

Ducharme-Barth
N.
,
Gruss
A.
,
Vincent
M.
,
Kiyoofuji
H.
,
Aoki
Y.
,
Pilling
G.
,
Hampton
J.
, et al.
2022
.
Impacts of fishery-dependent spatial sampling patterns on catch-per-unit-effort standardization: a simulation study and fishery application
.
Fisheries Research
,
246
:
106169
.

FAO
.
2021
.
FAO Yearbook. Fishery and Aquaculture Statistics 2019
.
FAO
,
Rome
.

Fenberg
P.
,
Roy
K.
2008
.
Ecological and evolutionary consequences of size-selective harvesting: how much do we know?
.
Molecular Ecology
,
17
:
209
220
.

Fenske
K.
,
Hulson
P.
,
Williams
B.
,
O'Leary
C
.
2020
.
Assessment of the Dusky Rockfish stock in the Gulf of Alaska
.
(last accessed 26 December 2022)
.

Francis
R.I.C.
2011
.
Data weighting in statistical fisheries stock assessment models
.
Canadian Journal of Fisheries and Aquatic Sciences
,
68
:
1124
1138
.

Gardner
J.
,
Peters
A.
,
Kearney
M.
,
Joseph
L.
,
Heinsohn
R.
2011
.
Declining body size: a third universal response to warming?
.
Trends in Ecology & Evolution
,
26
:
285
291
.

Gruss
A.
,
Thorson
J.T.
,
Stawitz
C.C.
,
Reum
J.C.P.
,
Rohan
S.K.
,
Barnes
C.L.
2021
.
Synthesis of interannual variability in spatial demographic processes supports the strong influence of cold-pool extent on eastern Bering Sea walleye pollock (Gadus chalcogrammus)
.
Progress in Oceanography
,
194
:
102569
.

Gruss
A.
,
Walter
J.
,
Babcock
E.
,
Forrestal
F.
,
Thorson
J.
,
Lauretta
M.
,
Schirripa
M.
2019
.
Evaluation of the impacts of different treatments of spatio-temporal variation in catch-per-unit-effort standardization models
.
Fisheries Research
,
213
:
75
93
.

Gunderson
D. R.
1993
.
Surveys of Fisheries Resources
.
John Wiley & Sons
,
New York, NY
.

Hanselman
D.
,
Heifetz
J.
,
Echave
K.
,
Dressel
S.
2015
.
Move it or lose it: movement and mortality of sablefish tagged in Alaska
.
Canadian Journal of Fisheries and Aquatic Sciences
,
72
:
238
251
.

Helser
T.
,
Brodziak
J.
1998
.
Impacts of density-dependent growth and maturation on assessment advice to rebuild depleted US siler hake (Merluccius bilinearis) stocks
.
Canadian Journal of Fisheries and Aquatic Sciences
,
55
:
882
892
.

Hilborn
R.
,
Walters
C. J.
1992
.
Quantitative Fisheries Stock Assessment: Choice, Dynamics and Uncertainty
, 1st edn.
Springer
,
Norwell, MA
.

Holsman
K.
,
Aydin
K.
,
Sullivan
J.
,
Hurst
T.
,
Kruse
G.
2019
.
Climate effects and bottom-up controls on growth and size-at-age of Pacific halibut (Hippoglossus stenolepis) in Alaska (USA)
.
Fisheries Oceanography
,
28
:
345
358
.

Hunt
G.
,
Stabeno
P.
,
Strom
S.
,
Napp
J.
2008
.
Patterns of spatial and temporal variation in the marine ecosystem of the southeastern Bering Sea, with special reference to the Pribilof Domain
.
Deep Sea Research Part II: Topical Studies in Oceanography
,
55
:
1919
1944
.

Hunt
G.
,
Yasumiishi
E.
,
Eisner
L.
,
Stabeno
P.
,
Decker
M.
2022
.
Climate warming and the loss of sea ice: the impact of sea-ice variability on the southeastern Bering Sea pelagic ecosystem
.
ICES Journal of Marine Science
,
79
:
937
953
.

Ianelli
J.
,
Fissel
B.
,
Stienessen
S.
,
Honkalehto
T.
,
Siddon
E.
,
Allen-Akselrud
C
.
2021
.
Chapter 1: Assessment of the Walleye Pollock Stock in the Eastern Bering Sea
.
(last accessed 26 December 2022)
.

Ikpewe
I.
,
Baudron
A.
,
Ponchon
A.
,
Fernandes
P.
2021
.
Bigger juveniles and smaller adults: changes in fish size correlated with warming seas
.
Journal of Applied Ecology
,
58
:
847
856
.

Jorgensen
C.
,
Ernande
B.
,
Oyvind
F.
2009
.
Size-selective fishing gear and life history evolution in the northeast Arctic cod
.
Evolutionary Applications
,
2
:
356
370
.

Kimura
D.
,
Somerton
D.
2006
.
Review of statistical aspects of survey sampling for marine fisheries
.
Reviews in Fisheries Science
,
14
:
245
283
.

Kotwicki
S.
,
Ianelli
J.
,
Punt
A.
2014
.
Correcting density-dependent effects in abundance estimates from bottom-trawl surveys
.
ICES Journal of Marine Science
,
71
:
1107
1116
.

Kotwicki
S.
,
Ono
K.
2018
.
The effect of random and density-dependent variation in sampling efficiency on variance of abundance estimates from fishery surveys
.
Fish and Fisheries
,
20
:
760
774
.

Kraak
S.
,
Haase
S.
,
Minto
C.
,
Santos
J.
2019
.
The Rosa Lee phenomenon and its consequences for fisheries advice on changes in fishing mortality or gear selectivity
.
ICES Journal of Marine Science
,
76
:
2179
2192
.

Krohn
M.
,
Kerr
S.
1997
.
Declining weight-at-age in northern cod and the potential importance of the early years and size-selective fishing mortality
.
Science Council Studies
,
29
:
43
50
.

Kuriyama
P.
,
Ono
K.
,
Hurtado-Ferro
F.
,
Hicks
A.
,
Taylor
I.
,
Licandeo
R.
,
Johnson
K.
, et al.
2016
.
An empirical weight-at-age approach reduces estimation bias compared to modeling parametric growth in integrated, statistical stock assessment models when growth is time varying
.
Fisheries Research
,
180
:
119
127
..

Lindgren
F.
,
Rue
H.
2015
.
Bayesian spatial modelling with R-INLA
.
Journal of Statistical Software
,
63
:
1
25
.

Lindgren
F.
,
Rue
H.
,
Lindstrom
J.
2011
.
An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach
.
Journal of the Royal Statistical Society: Series B (Statistical Methodology)
,
73
:
423
498
.

Lindmark
M
,
Ohlberger
J.
,
Gardmark
A.
2022
.
Optimum growth temperature declines with body size within fish species
.
Global Change Biology
,
28
:
2259
2271
.

Link
J.
,
Nye
J.
,
Hare
J.
2011
.
Guidelines for incorporating fish distribution shifts into a fisheries management context
.
Fish and Fisheries
,
12
:
461
469
.

Martin
T.
,
Wintle
B.
,
Rhodes
J.
,
Kunert
P.
,
Field
Sa.
,
Low-Choy
S.
,
Tyre
A.
, et al.
2005
.
Zero tolerance ecology: improving ecological inference by modelling the source of zero observations
.
Ecology Letters
,
8
:
1235
1246
.

Maunder
M.
,
Punt
A.
2004
.
Standardizing catch and effort data: a review of recent approaches
.
Fisheries Research
,
70
:
141
159
.

Maunder
M.
,
Thorson
J.
,
Xu
H.
,
Oliveros-Ramos
R.
,
Hoyle
S.
,
Tremblay-Boyer
L.
,
Lee
H.
, et al.
2020
.
The need for spatio-temporal modeling to determine catch-per-unit effort based indices of abundance and associated composition data for inclusion in stock assessment models
.
Fisheries Research
,
229
:
105594
.

Methot
R.
2009
.
Stock assessment: operational models in support of fisheries management
. In
The Future of Fisheries Science in North America
, pp.
137
165
.. Ed. by
R.J.
Beamish
,
B.J.
Rothschild
.
Springer
,
Dordrecht
.

Minte-Vera
C.
2004
.
Meta-analysis of density-dependent growth
.
PhD thesis
,
University of Washington
,
Seattle, WA
,
1
204
.pp.

Mize
T.
2019
.
Best practices for estimation, interpreting, and presenting nonlinear interaction effects
.
Sociological Science
,
6
:
81
117
.

Mueter
F. J.
,
Bond
N.A.
,
Ianelli
J.N.
,
Hollowed
A.B.
2011
.
Expected declines in recruitment of walleye pollock (Theragra chalcogramma) in the eastern Bering Sea under future climate change
.
ICES Journal of Marine Science
,
68
:
1284
1296
.

Mueter
F.
,
Litzow
M.
2008
.
Sea ice retreat alters the biogeography of the Bering Sea continental Shelf
.
Ecological Applications
,
18
:
309
320
.

National Oceanic and Atmospheric Administration
.
2020a
.
NOAA Fisheries Will Cancel Five Alaska Research Surveys for 2020
.
Media Release
.
(last accessed 26 December 2022)
.

National Oceanic and Atmospheric Administration
.
2020b
.
NOAA Fisheries Cancels 2020 Southeast Fishery-Independent Survey (SEFIS)
.
Media Release
.
(last accessed 26 December 2022)
.

National Oceanic and Atmospheric Administration
.
2020c
.
NOAA Fisheries Cancels Three West Coast Surveys for 2020
.
Agency Statement
.
(last accessed 26 December 2022)
.

National Oceanic and Atmospheric Administration
.
2020d
.
NOAA Fisheries Cancels Three West Coast Surveys for 2020
.
Agency Statement
.
(last accessed 26 December 2022)
.

Nebenzahl
D.
,
Goddard
P.
(Compilers)
.
2000
.
Bottom trawl survey of the eastern Bering Sea continental shelf
.
AFSC Processed Rep. 2000-10. Alaska Fish. Sci. Cent., NOAA, Natl. Mar, Fish. Serv., 7600 Sand Point Way NE, Seattle WA
.

O'Leary
C.
,
DeFilippo
L.
,
Thorson
J.
,
Kotwicki
S.
,
Hoff
G.
,
Kulik
V.
,
Ianelli
J.
, et al.
2022
.
Understanding transboundary stocks’ availability by combining multiple fisheries-independent surveys and oceanographic conditions in spatiotemporal models
.
ICES Journal of Marine Science
,
79
:
1063
1074
.

O'Leary
C.
,
Thorson
J.
,
Ianelli
J.
,
Kotwicki
S.
2020
.
Adapting to climate-driven distribution shifts using model-based indices and age composition from multiple surveys in the walleye pollock (Gadus chalcogrammus) stock assessment
.
Fisheries Oceanography
,
29
:
541
557
.

Oke
K.
,
Mueter
F.
,
Litzow
M.
2022
.
Warming leads to opposite patterns in weight-at-age for young versus old age classes of Bering Sea walleye pollock
.
Canadian Journal of Fisheries and Aquatic Sciences
,
79
:
1655
1666
.,

Pinsky
M.
,
Selden
R.
,
Kitchel
Z.
2020
.
Climate-driven shifts in marine species ranges: scaling from organisms to communities
.
Annual Review of Marine Science
,
12
:
153
179
.

Punt
A.
,
Smith
D.
2001
.
Assessments of species in the Australian South East Fishery can be sensitive to the method used to convert from size- to age-composition data
.
Marine & Freshwater Research
,
52
:
683
690
.

Punt
A.E.
,
Dalton
M. G.
,
Cheng
W.
,
Hermann
A.J.
,
Holsman
K.K.
,
Hurst
T.P.
,
Ianelli
J. N.
, et al.
2021
.
Evaluating the impact of climate and demographic variation on future prospects for fish stocks: an application for northern rock sole in Alaska
.
Deep Sea Research Part II: Topical Studies in Oceanography
,
189–190
:
104951
.

Rooper
C.
,
Ortiz
I.
,
Hermann
A.
,
Laman
N.
,
Cheng
W.
,
Kearney
K.
,
Aydin
K.
2021
.
Predicted shifts of groundfish distribution in the eastern Bering Sea under climate change, with implications for fish populations and fisheries management
.
ICES Journal of Marine Science
,
78
:
220
234
.

Sharpe
D.
,
Hendry
A.
2009
.
SYNTHESIS: Life history change in commercially exploited fish stocks: an analysis of trends across studies
.
Evolutionary Applications
.
2
(
3
):
260
275
.. https://doi.org/10.1111/j.1752-4571.2009.00080.x.

Sullivan
J.
2016
.
Environmental, ecological, and fishery effects on growth and size-at-age of Pacific halibut (Hippoglossus stenolepis)
.
MSc thesis
,
University of Alaska Fairbanks
,
Fairbanks, AK
.

Sullivan
J.
,
Kruse
G.
,
Mueter
F.
2018
.
Do environmental and ecological conditions explain declines in size-at-age of Pacific halibut in the Gulf of Alaska
. In
Impacts of a Changing Environment on the Dynamics of High-latitude Fish and Fisheries
. Ed. by
F.J.
Mueter
,
M.R.
Baker
,
S.C.
Dressel
,
A.B.
Hollowed
.
Alaska Sea Grant, University of Alaska Fairbanks
,
Fairbanks, AK
.

Swain
D.
,
Sinclair
A.
,
Hanson
J.
2007
.
Evolutionary response to size-selective mortality in exploited fish population
.
Proceedings of the Royal Society B: Biological Sciences
,
274
:
1015
1022
.

Thorson
J.
2014
.
Standardizing compositional data for stock assessment
.
ICES Journal of Marine Science
,
71
:
1117
1128
.

Thorson
J.
2018
.
Three problems with the conventional delta-model for biomass sampling data, and a computationally efficient alternative
.
Canadian Journal of Fisheries and Aquatic Sciences
,
75
:
1369
1382
.

Thorson
J.
2019
.
Guidance for decisions using the Vector Autoregressive Spatio-Temporal (VAST) package in stock, ecosystem, habitat and climate assessments
.
Fisheries Research
,
210
:
143
161
.

Thorson
J.
,
Haltuch
M.
2019
.
Spatiotemporal analysis of compositional data: increased precision and improved workflow using model-based inputs to stock assessment
.
Canadian Journal of Fisheries and Aquatic Sciences
,
76
:
401
414
.

Thorson
J.
,
Hermann
A.
,
Siwicke
K.
,
Zimmermann
M.
2021
.
Grand challenge for habitat science: stage-structured responses, nonlocal drivers, and mechanistic associations among habitat variables affecting fishery productivity
.
ICES Journal of Marine Science
,
78
:
1956
1968
.

Thorson
J.
,
Ianelli
J. N.
,
Kotwicki
S.
2017
.
The relative influence of temperature and size-structure on fish distribution shifts: a case-study on walleye pollock in the Bering Sea
.
Fish and Fisheries
,
18
:
1073
1084
.

Thorson
J.
,
Maunder
M.
,
Punt
A.
2020
.
The development of spatio-temporal models of fishery catch-per-unit-effort data to derive indices of relative abundance
.
Fisheries Research
,
230
:
105611
.

Thorson
J.
,
Minte-Vera
C.
2016
.
Relative magnitude of cohort, age, and year effects on size at age of exploited marine fishes
.
Fisheries Research
,
180
:
45
53
.

Thorson
J.
,
Monnahan
C.
,
Cope
J.
2015a
.
The potential impact of time-variation in vital rates on fisheries management targets for marine fishes
.
Fisheries Research
,
169
:
8
17
.

Thorson
J.
,
Shelton
A.O.
,
Ward
E.J.
,
Skaug
H.J.
2015b
.
Geostatistical delta-generalized linear mixed models improve precision for estimated abundance indices for West Coast groundfishes
.
ICES Journal of Marine Science
,
72
:
1297
1310
.

Thorson
J.
,
Stewart
I.
,
Punt
A.
2011
.
Accounting for fish shoals in single- and multi-species survey data using mixture distribution models
.
Canadian Journal of Fisheries and Aquatic Sciences
,
68
:
1681
1693
.

Uusi-Heikkila
S.
,
Whiteley
A.
,
Kuparinen
A.
,
Matsumura
S.
,
Venturelli
P.
,
Wolter
C.
,
Slate
J.
, et al.
2015
.
The evolutionary legacy of size-selective harvesting extends from genes to populations
.
Evolutionary Applications
,
8
:
597
620
.

van Rijn
I.
,
Buba
Y.
,
DeLong
J.
,
Kiflawi
M.
,
Belmaker
J.
2017
.
Large but uneven reduction in fish size across species in relation to changing sea temperatures
.
Global Change Biology
,
23
:
3667
3674
.

Walters
G.
1997
.
1993 Bottom Trawl Survey of the Eastern Bering Sea Continental Shelf. AFSC Processed Rep. 97-09. Alaska Fish. Sci. Cent., NOAA, Natl. Mar, Fish. Serv., 7600 Sand Point Way NE, Seattle WA
.

Whitehouse
G.
,
Aydin
K.
,
Hollowed
A.
,
Holsman
K.
,
Cheng
K.
,
Faig
A.
,
Haynie
A.
, et al.
2021
.
Bottom-up impacts of forecasted climate change on the eastern Bering Sea food web
.
Frontiers in Marine Science
,
8
:
624301
.

Whitten
A.
,
Klaer
N.
,
Tuck
G.
,
Day
R.
2013
.
Accounting for cohort-specific variable growth in fisheries stock assessments: a case study from south-eastern Australia
.
Fisheries Research
,
142
:
27
36
.

Williams
A.
,
Davies
C.
,
Russ
G.
2003
.
Scales of spatial variation in demography of a large coral-reef fish-an exception to the typical model?
.
Fisheries Bulletin
.
101
:
673
683
.

Williams
R.
2012
.
Using the margins command to estimate and interpret adjusted predictions and marginal effects
.
The Stata Journal
,
12
:
308
331
.

Wilson
M.
,
Armistead
C.
1991
.
NOAA Technical Memorandum NMFS F/NWC-212
.
Alaska Fish. Sci. Cent., NOAA, Natl. Mar, Fish. Serv., 7600 Sand Point Way NE, Seattle WA
.

Zhou
S.
,
Campbell
R.
,
Hoyle
S.
2019
.
Catch per unit effort standardization using spatio-temporal models for Australia’s eastern Tuna and Billfish Fishery
.
ICES Journal of Marine Science
,
76
:
1489
1504
.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
Handling Editor: Pamela Woods
Pamela Woods
Handling Editor
Search for other works by this author on:

Supplementary data