Assessing benthos through predator stomach contents: spatiotemporal modeling of abundance and habitat use

Abstract

In shelf ecosystems, benthic invertebrates facilitate nutrient recycling and the transfer of energy to higher trophic levels. However, large‐scale monitoring through direct sampling (e.g. using benthic grabs or bottom trawls) can be costly in terms of time and labor. Here, we demonstrate a method for developing standardized abundance indices of forage groups (i.e. species or functional groups preyed upon by predators) based on predator stomach contents. The modeling approach is analogous to methods for estimating abundance indices from fisheries catch‐per‐unit‐effort data; accounts for predator species‐ and size‐specific differences in forage group selectivities, which may vary over space; and permits index estimation when diets are unevenly sampled over space, time, and across predator species. We apply the model to four decades of groundfish diet data from the eastern Bering Sea and estimate abundance indices for nine benthic forage groups. The fitted models were then used to generate: 1) time‐averaged maps of relative biomass density for forage groups and demarcation of potential core habitat areas; 2) region‐wide biomass index time series; and 3) an assessment of bioregions based on benthic forage community composition and change in bioregion area over time. Diet‐based biomass densities were on average correlated well with densities obtained from direct sampling (bottom trawl) for species of Chionoecetes crabs (0.61–0.69). Correlations for polychaetes and bivalves, which had fewer direct survey samples (benthic grabs) were also positive but weaker (0.45 and 0.20, respectively), potentially reflecting larger differences in selectivities between predators and sampling gear or sampling error. We argue that diet data can provide an additional and cost‐efficient lens through which abundances of forage species can be quantified and aid efforts to monitor change in marine ecosystems. Abundance indices can be used in subsequent whole‐of‐ecosystem models, and habitat utilization maps can be used to inform spatial management decisions.

Introduction

Benthic invertebrates are a major component of marine food webs, facilitate nutrient recycling, and provide critical forage for fished species (Grebmeier et al. 1988, Reynoldson and Metcalfe-Smith 1992, Hale et al. 2017). For example, in the eastern Bering Sea (EBS), a highly productive subpolar shelf ecosystem, they account for 58% of animal biomass and support 52% of fish production (Aydin et al. 2007, Whitehouse et al. 2021). Accordingly, they may act as ecological indicators, and changes in their abundances may help reveal shifts in food web structure, function, and productivity (Grebmeier et al. 2006, Yeung and McConnaughey 2006). However, monitoring benthic invertebrate abundances in marine ecosystems through direct sampling (e.g. benthic grabs, bottom trawls, video surveys) can be expensive in terms of time and labor (Bilyard 1987, Carey and Keough 2002), constraining the spatial and temporal scales at which sampling is conducted. Moreover, different sampling gears are effective at capturing different subsets of the benthic community (Jørgensen et al. 2011), and use of only one gear may miss forage groups and size classes fed upon by higher trophic level functional groups.

A promising, complementary approach to sampling with mechanical gear entails the treatment of predators as biological samplers (Link 2004). Predator stomach content data provide a ‘snapshot' of recent foraging success and the prevalence of forage groups in diets can indicate availability and abundance (Link 2004, Yeung et al. 2010, Deroba 2018, Ng et al. 2021, Gaichas et al. 2024). If predators are routinely sampled as part of ongoing resource monitoring efforts, the collection of diet data is relatively straightforward and the resulting forage abundance index adopts the spatial and temporal scale of the predator distribution within the resource survey (Ng et al. 2021). Further, by their very nature, diet-based indices directly characterize abundances of forage groups and sizes that are trophically relevant and therefore potentially important indicators of ecosystem change. That said, inferring forage abundance trends from diet data presents challenges. For instance, prey taxonomy may be difficult to resolve depending on digestion level, and more challenging for soft-bodied prey (e.g. polychaetes) which are easier to digest relative to hard-bodied prey (e.g. crustaceans; Andersen et al. 2016). If characterizing abundances over large areas, the synthesis of diet data from multiple predator species may also be required, thus necessitating methods that account for predator- and size-specific preferences for forage species or functional groups (Deroba 2018). Within predator species, forage preference or catchability may also change as a function of space due to habitat features or the availability of alternative prey (Ng et al. 2021, Grüss et al. 2023).

Recent research describes stomach content data as if they arise from the predator capturing some proportion of prey that it encounters within a given foraging area, where individual prey are spatially distributed according to their habitat preferences and have unique body sizes. These assumptions result in a thinned and double-marked point process (Thorson et al. 2022). In more detail, this process arises when individual prey have a random location that follows some underlying density function (i.e. following a point process), where each individual prey has some body size (a continuous mark), and a taxonomic or functional label (a categorical mark). Next, envision that the contents of each stomach correspond to a foraging trip over some portion of this landscape with unknown area, and that diet samples arise from predator foraging (i.e. attack, capture, handling, and digestion) rates over this foraging area. Under this model, diet samples can be approximated using a multivariate generalized linear mixed model (GLMM), where the expected quantity of each forage category arises from the product of forage numerical density, body size, and predator foraging rates. Expressing stomach contents as a GLMM allows us to estimate how the product of numerical density, average size, and foraging rates varies across space and time. Foraging rates then are interpreted similarly to ‘gear catchability' (a.k.a. detectability) in conventional sampling theory, and are either assumed to be constant over space (when estimating habitat utilization) or time (when estimating abundance indices), and can also be intercalibrated or validated relative to other sampling gears.

Building on this conceptual model for foraging processes and their treatment as a GLMM, spatiotemporal models have emerged as promising tools for synthesizing and standardizing ecological data sets including predator stomach content data (Grüss et al. 2021, 2023, Ng et al. 2021, Gaichas et al. 2024). For instance, using diet data collected from six different Northwest Atlantic predator species, Ng et al. (2021) developed abundance indices for the forage species Atlantic herring Clupea harengus wherein separate spatiotemporal models were fit for each predator in which the specific mass of herring (that is, g prey g⁻¹ predator) was treated as the response variable. Separate abundance indices based on each predator were obtained by summing the specific mass of Atlantic herring (controlling for predator size) across the survey region within years. To better synthesize data from different sources into a common index, Grüss et al. (2023) adopted an approach that differed in two significant ways. First, data from multiple sources were included in a single model (in their case, they included snow crab Chionoecetes opilio data from different surveys as well as diet data from Pacific cod Gadus macrocephalus, a predator of snow crab) along with covariates to account for differences in catchability between data sources. Second, covariates specific to each data source were allowed to vary spatially to account for differences in prey catchability. The resulting abundance index incorporated information from all data sources while better accounting for systematic differences in prey mass or catches. The method is analogous to those used to construct abundance indices based on catch-per-unit-effort (CPUE) when different surveys or sampling gears are combined (Maunder and Punt 2004, Grüss et al. 2023). The extensions proposed by Grüss et al. (2023) offer a potentially useful way to combine diet data from multiple predators to generate abundance indices of forage species, but they are relatively new and have seen limited application.

The EBS is a large, subpolar ecosystem with seasonal ice coverage that overlies a broad continental shelf (Fig. 1A). During the spring and summer growing season, high rates of primary production result in large fluxes of detritus to the seafloor which support high benthic standing stocks and productivity levels (Stoker 1978, Grebmeier et al. 1988, 2015, Aydin et al. 2007, Whitehouse et al. 2021). In the adjacent northern Bering Sea (NBS), research programs have undertaken long-term benthic monitoring at locations in the vicinity of St. Lawrence Island (Grebmeier 2015, Moore and Grebmeier 2018). Similarly, in the EBS, summertime surveys since the 1970s have provided important information on spatial patterns and temporal trends in groundfish and large-bodied epibenthic invertebrates vulnerable to bottom trawls (Stoker 1978, Jewett and Feder 1981, Yeung and McConnaughey 2006, Stevenson and Lauth 2012). However, within the EBS little attention has been paid to smaller epibenthos and infauna which comprise larger percentages of benthic invertebrate biomass (30 and 61 versus 9%, respectively; Whitehouse et al. 2021). Large-scale surveys (i.e. covering areas > 1000 km²) that directly sampled for infauna and epibenthos are limited to the late 1950s, 1970s, and mid-1980s (Haflinger 1981, McDonald et al. 1981, Stoker 1981, Grebmeier et al. 1988, Coyle et al. 2007). Since 2006, smaller-scale surveys have also occurred on a sporadic basis, but they are concentrated in the southern EBS and have limited replication over time and space due to the changing nature of the research objectives motivating the sampling (Yeung et al. 2010, Yeung and Yang 2017, 2017). Collectively, these studies have resolved important infaunal and epibenthos patterns including overall increasing trends in biomass with latitude and relationships between community composition and covariates such as distance to shore, depth, substrate characteristics, and other spatially structured variables including properties of overlaying water masses. The sparsity of sampling with respect to time, however, constrains efforts to characterize temporal variation and system responses to environmental conditions. In contrast, demersal fishes in the region have been sampled using a standardized gear and protocol nearly continuously since 1982 in an annual, region-wide fishery-independent bottom trawl survey (Stevenson and Lauth 2012). Although the gear is ill-suited for directly sampling small-bodied epibenthic invertebrates (i.e. less than 30 mm body width; Kotwicki et al. 2017), food habits data are collected from several predator species, many of which feed heavily on benthic forage groups including infauna and epibenthos (Lang et al. 2003, Livingston et al. 2017).

[IMAGE OMITTED. SEE PDF]

Here, we sought to synthesize diet data from multiple predators using spatiotemporal models to develop region-wide abundance indices for key EBS benthic forage groups. In doing so, we aimed to maximize the spatial and temporal extent of data coverage contributing to each index. Specifically, we estimated spatiotemporal patterns of benthic forage groups that correspond to functional groups or nodes in EBS food web models (Aydin et al. 2007, Whitehouse et al. 2021) and are therefore of interest with respect to modeling trophic flows. Annual abundance indices of the forage groups may thus also hold value for characterizing ecosystem trends and contribute to ecosystem indicator development efforts (Zador et al. 2017, Siddon 2024). In addition, we demonstrate how the fitted models can be used to generate a wider range of products with relevance for spatial management and the characterization of biogeographic patterns. Although we focus here on the estimation of spatiotemporal patterns, we note that the modeling framework can also be used to test causal hypotheses regarding forage distributions, and is the subject of ongoing research.

Material and methods

Our analysis proceeded in two phases. First, we fit spatiotemporal models to the diets of benthic predators to estimate densities for several important benthic forage groups. Where possible, we evaluated the level of correlation between diet-based forage density estimates and those obtained from direct survey-based measurements. In the second phase, we used the fitted models to estimate 1) distribution maps and assess core benthic forage group habitat areas, 2) obtain region-wide annual biomass indices, and 3) evaluate bioregions in the EBS based on benthic forage community composition and change in bioregion area over time. We first provide overviews of the EBS diet data and modeling framework, then describe analytical approaches for achieving objectives in the first and second analysis phases.

Eastern Bering Sea diet data

We obtained diet data for 11 predator species that feed significantly on benthic forage groups (Supporting information). All were sampled as part of a fishery-independent bottom trawl summer groundfish survey of the EBS by NOAA's Alaska Fisheries Science Center (Livingston et al. 2017). The survey covers the southeastern and northern Bering Sea (SEBS and NBS), with areas of 4.9 × 10⁵ and 2.0 × 10⁵ km², respectively. Sampling stations are gridded, with 20 km spacing (Fig. 1A; Stevenson and Lauth 2012). In the SEBS, the survey has been performed nearly continuously using the same sampling design since 1982. For the NBS, the first survey was performed in 2010, again in 2017 and 2019, and annually from 2021 to 2023. As part of the surveys, stomach content data are collected from predators using an area- and length-stratified sampling design (Livingston et al. 2017). We focused our analyses on diet data collected between 1984 and 2022, when sampling intensity was generally higher and included multiple predators, with the exception of years when the bottom trawl survey was not performed (2020) or sufficient diet data were not collected (2004 and 2020; Fig. 1B). Additional details regarding the bottom trawl survey design and stomach sampling are provided in the Supporting information. All analyses were limited to 1984 through 2022 based on diet data availability (Fig. 1B).

Based on previous food web modeling efforts (Aydin et al. 2007, Whitehouse et al. 2021), we identified nine benthic invertebrate forage groups that, combined, support 45% of consumption demands for all EBS fish: polychaetes, benthic amphipods, benthic shrimps, bivalves, brittle stars, snow crab C. opilio, Tanner crab Chionoecetes bairdi, miscellaneous worms, and miscellaneous crustaceans. The group ‘miscellaneous worms' primarily consists of annelid worms (including oligochaetes, leeches, flatworms) and sipunculids while ‘miscellaneous crustaceans' consists of epifauna such as cumaceans, isopods, and small-bodied arthropods (see the Supporting information for additional taxonomic details of each forage group). Although most of the forage groups are taxonomically diverse, they correspond to distinct functional groups or nodes in regional food web models (Aydin et al. 2007, Whitehouse et al. 2021) and are therefore relevant for system monitoring in the EBS.

Predators were included as samplers of a forage group if 1) diet records were available from at least 1000 individuals and 2) the forage group appeared in at least 10% of individual diets. The number of predators included in forage group models ranged between seven (polychaetes and benthic amphipods) and one (snow and Tanner crabs; Supporting information). Maps depicting stomach sampling intensity over the survey area are provided for each predator species in the Supporting information. For two of the most widely distributed predators in the region, Pacific cod and walleye pollock, the proportion of the survey area covered by stomach samples ranged between ~ 80 and 97% across most years (Supporting information). For the remaining predators, spatial coverage was lower and more variable due to changes in sampling effort over time and differences in predator distributions (Supporting information). However, when stomach samples from all potential samplers of a forage group were pooled, spatial coverage of the survey area was higher overall and typically greater than 80% for all but three groups (brittle stars, bivalves, and miscellaneous crustaceans). For those groups, spatial coverage was relatively high (> 70%) in years prior to 2000, but declined markedly thereafter (Supporting information).

For each individual predator, we first obtained the specific forage mass (g forage g⁻¹ predator), where predator mass was calculated from predator length using published length–weight relationships (Reum et al. 2019). Predators were then binned according to 5 cm length intervals, and specific forage masses were averaged by sample date, station, and size bin. The resulting mean value was treated as the response variable in statistical models. Averaging was deemed necessary to reduce the data set size and improve computational efficiency of the analysis (Grüss et al. 2021). Predator length, which is included as a covariate in models, was obtained in a similar manner.

Spatiotemporal models

We sought to estimate prey biomass densities as a thinned and double-marked point process. To do so, we fit the model using the multivariate spatiotemporal package ‘tinyVAST' (Thorson et al. 2024a). Here, we describe the model using vector-matrix notation, i.e. bold for vectors (Edwards and Auger-Méthé 2019), but modify notation from ‘tinyVAST' for clarity of presentation.

We specify a log-linked linear predictor for prey biomass density d_s,t at discretized locations s within a continuous spatial domain $s\in \mathcal{D}$ and annual time intervals $t\in \left\{1982,1983\dots 2022\right\}$: $\log \left(d_{s,t}\right)={\beta }_t+{\omega }_s+{\epsilon }_{s,t}$

where β_t is an annually varying intercept, ω is a spatially correlated random effect that is constant over time: $$

and where R is a spatial correlation matrix that approximates a Matérn correlation function that accounts for geometric anisotropy (Thorson et al. 2015) and uses the SPDE method (Lindgren et al. 2011) with estimated decorrelation rate κ, ${\sigma }_{\omega }^2$ is the estimated magnitude of spatial variance, and MVN indicates a multivariate normally distributed variable. In practice, ω (and other spatially correlated random effects in the model) are represented using a triangulated mesh, where vertices (or ‘knots') are used to approximate spatial variability of the observations and bilinear interpolation is used to generate a continuous spatial field.

Similarly, ε_t is a spatially correlated random effect that varies for each year t: $$

where ${\sigma }_{\epsilon }^2$ is the estimated spatiotemporal variance and ρ is the estimated magnitude of first-order autoregression in the spatiotemporal term (where ω and ε_t use the same estimated decorrelation rate κ).

However, we estimate prey biomass density d_s,t using biomass samples y_i from C predators, where c[i] indicates the predator for a given sample i. We assume that biomass consumption can be approximated as foraging individual prey items (i.e. approximately a Poisson point process) where each prey item follows a continuous size (i.e. approximately a gamma distribution). These assumptions result in a Tweedie distribution with intensity λ_i for each sample (Thorson et al. 2022), which can be fitted as a generalized linear mixed model: $y_i\sim \text{Tweedie}\left({\lambda }_i,{\varphi }_{c[i]},p\right)$

where ϕ_c[i] is the dispersion for each predator c and p is the power parameter that is shared across predators. Similarly, foraging intensity λ_i is affected by prey biomass density as well as an additional term representing thinning rates, i.e. variation in foraging area as well as attack, capture, and digestion rates among predators: $$

where catchability q_i includes the expected log-ratio γ_c[i] for foraging rates for each predator c, as well as a spatially varying coefficient (SVC) ξ_s[i],c[i] (Thorson et al. 2023) representing variation in foraging rates for each predator c at location s: $\log \left(q_i\right)=\underset{\text{intercept}}{\underbrace{{\gamma }_{c[i]}}}+\underset{\text{SVC}}{\underbrace{{\xi }_{s[i],c[i]}}}+\underset{\text{covariates}}{\underbrace{{\alpha }_{1,c[i]}{x}_{s,t}+{\alpha }_{2,c[i]}{x}_{s,t}^2+{\delta }_{c[i]}{l}_i}}$

and where SVC ξ_s,c for predator c at location s $$

where ${\tau }_c^2$ is the estimated SVC variance for each predator. For identifiability, we specify that ${\gamma }_{c^{*}}=0$ and ${\xi }_{s,c^{*}}=0$ for a selected reference predator c^*, such that γ_c for other predators is their expected ratio relative to reference predator c^*. We also estimate a dome-shaped response of foraging rates to bottom temperatures for each predator by including a separate response (α_1,c and α_2,c) for each predator to bottom temperature x_s,t and bottom-temperature squared $x_{s,t}^2$ as covariates. Finally, we include a separate response δ_c for each predator to predator body length l_i for each diet sample as covariate to account for size-based shifts in foraging rates and preference. We center temperature and predator length (by subtracting off their mean) prior to model fitting.

For all forage groups, the model domain spanned both the SEBS and NBS, and spatial components of the model used a mesh with a minimum distance between knots (i.e. cutoff) of 60 km, resulting in 352 knots total. To fit the models, ‘tinyVAST' utilizes Template Model Builder, which employs the Laplace method to approximate the marginal likelihood. It then optimizes the log-marginal likelihood using gradients computed using automatic differentiation (Kristensen et al. 2016). After fitting, model convergence was confirmed (all marginal log-likelihood gradients were < 1.0 × 10⁻⁴ for fixed effects and the Hessian matrix of secondary derivatives of the negative likelihood were positive-definite). To evaluate model fit, we used simulation-based probability-integral-transform (PIT) residuals (Smith 1985, Warton et al. 2017) visualized using R package ‘DHARMa' (Hartig and Lohse 2020) and examined histograms and QQ-plots of the ‘DHARMa' residuals. All forage group models successfully converged and exhibited normally distributed residuals and relatively uniform distributions of PIT residuals between 0 and 1, indicative of adequate model fit (Hartig and Lohse 2020, Supporting information).

Composite measures of local prey densities

Under the modeling framework, the SVC term ξ_s[i],c[i] permits foraging rates to differ over space between predator species. Therefore, maps of relative forage densities depend on the SVC effect sizes of the specified predator. If a rationale exists for selecting one predator species as a better sampler of a forage group over another, distribution maps based on that predator alone may be appropriate (Grüss et al. 2023). However, analysts may lack a strong rationale for doing so. Instead, composite measures of forage density that combine predictions from all predators may be appropriate. We considered two possible approaches. For the first, we assumed that predator-specific spatial processes that manifest as different SVCs have an equal chance of increasing or decreasing predicted foraging rates relative to other predators in the model (for instance, the availability of alternative forage groups or presence of habitat features that change foraging success). We then defined the composite measure, D_Eq,s,t, as an average of the expected relative foraging rates of the predators, where: $D_{\text{Eq},s,t}=\frac{1}{C}{\displaystyle \sum }_{c=1}^C{\widehat{\lambda }}_{s,t,c}$

For each predator c we obtained expected foraging rates ${\widehat{\lambda }}_{s,t,c}$ for locations s in year t, and body lengths for each predator species that corresponded to the size that on average fed most heavily on the forage group. Specifically, for each predator we used their mean body length weighted by the specific mass of the forage group from the diet data set.

In the second composite measure, we made the additional assumption that the average prevalence of a forage group in the stomachs of a predator was proportional to its overall effectiveness as a sampler. That is, predators that fed heavily on a forage group were in general better samplers of that forage group. For the resulting composite measure, we obtained ${\widehat{\lambda }}_{s,t,c}$ as before, but weighted values for each predator by their mean expected foraging rate. The weighted composite measure, D_Wt,s,t, is given by: $D_{\text{Wt},s,t}=\frac{1}{C}{\displaystyle \sum }_{c=1}^C{\widehat{\lambda }}_{s,t,c}{w}_c$

where the normalized weights, w_c, for each predator c are calculated as: $w_c=\frac{{\widehat{\lambda }}_{\text{Ave},c}}{{\displaystyle {\sum }_{c=1}^C{\widehat{\lambda }}_{\text{Ave},c}}}$

and ${\widehat{\lambda }}_{\text{Ave},c}$ is the geometric mean of all ${\widehat{\lambda }}_{s,t,c}$ for predator c. We subsequently compared relative density maps based on D_Eq,s,t and D_Wt,s,t and a third measure based on expected foraging rates for the references predator alone (that is, simply ${\lambda }_{s,t,c}^{*}$ or D_R,s,t; Grüss et al. 2023).

Correspondence between diet- and survey-based density estimates

We evaluated diet-based estimates of relative forage densities against observed densities for four benthic forage groups with available data: snow and Tanner crabs, polychaetes, and bivalves. For crabs, we obtained biomass density estimates directly from the annual bottom trawl survey and calculated densities for crabs with carapace widths (CW) < 90 mm. The bottom trawl utilized in the survey has a 32-mm mesh codend liner and is considered useful for sampling Chionocetes crabs with CW > 30 mm (Kotwicki et al. 2017). Because Pacific cod, their main groundfish predator, feed on individuals between 20 and 90 mm CW, trawl density estimates were deemed suitable for validating diet-based estimates. For these forage species, the number of stations with both survey- and diet-based density estimates within a year was ~ 350. Due to the prevalence of zeros in the survey biomass density observations and failure to meet normality assumptions, we computed the spatial correlation within each year using the nonparametric Spearman's ranked correlation coefficient (called ‘spatial correlations' hereafter), and visualized the distribution of spatial correlations across all years. Although ranked correlations are less sensitive to outliers and linearity assumptions than parametric alternatives, they are best for characterizing monotonic relationships between variables. Trawl density estimates were not used to validate diet-based estimates for other epibenthic forage groups because the trawl is less effective at sampling organisms with significantly smaller body widths (e.g. benthic amphipods and shrimp) and those that are fragile and easily crushed by other catch in the net (e.g. brittle stars).

For polychaetes and bivalves, we obtained field density estimates based on samples collected using a 0.1 m² Van Veen benthic grab (Yeung et al. 2010, Yeung and Yang 2014, 2017, 2025). Benthic samples were available for 11 years between 2006 and 2019, but covered a substantially smaller proportion of the survey area in a given year (Supporting information). On average, 23 stations were sampled within a year, but ranged from 10 to 42; temporal replication across stations was low (Supporting information). We paired benthic grab biomass densities with diet-based composite density estimates averaged within a 10 km radius of the benthic grab location and calculated spatial correlations in the same manner as that for crabs.

Index standardization

Given the fitted model, we estimated area-expanded biomass indices for each forage group for the SEBS where annual sampling has occurred the longest. To do so, we divide the EBS model domain $\mathcal{D}$ into a set of 1248 square grid cells (each with an area of 10.4 km²), and record the spatial area a_g and centroid for each cell. We then predict expected foraging rates ${\widehat{\lambda }}_{s,t,c^{*}}$ at locations s in year t for the ‘reference predator' c^* at the sample-average predator length: $\log \left({\widehat{\lambda }}_{s,t,c^{*}}\right)={\widehat{\beta }}_t+{\widehat{\omega }}_s+{\widehat{\epsilon }}_{s,t}+{\widehat{\alpha }}_{1,c^{*}}{x}_{s,t}+{\widehat{\alpha }}_{2,c^{*}}{x}_{s,t}^2$

Next, we calculate the forage biomass index b_t as the area-expanded expected foraging for that predator: $b_t={\displaystyle \sum }_{s=1}^{n_s}{a}_s{\widehat{\lambda }}_{s,t,c^{*}}$

where a_s is the area associated with each prediction s. To mitigate bias arising from retransformation and skewness of random effects, we apply the epsilon estimator (Thorson and Kristensen 2016), and we also calculate the standard error using a generalization of the delta method (Kass and Steffey 1989). This model-based estimator is essentially ‘filtering out' the variation in predator foraging rates q_s,c by converting these to be in units of expected forage biomass for the reference predator c^* and at standard length. Due to the additive structure of the model and common spatiotemporal effect ε_t, the resulting index time series is often insensitive to the choice of predator as reference when expressed in anomaly space.

Forage distribution and hotspot maps

For each forage group, we generated time-averaged distribution maps of relative biomass densities for the entire EBS survey area (SEBS and NBS) based on composite density measures and evaluated evidence for forage biomass hotspots in terms of agreement across predators. For the latter, we defined biomass hotspots as the top fifth of EBS area with the highest relative forage densities. Given a forage group, we obtained time-averaged relative density maps for each predator species. Next, we categorized 20% of EBS area that accounted for the highest forage densities as a forage density hotspot for each predator. We noted areas where agreement across predator species in hotspot designation was greater than 50 and 80%.

Bioregion analysis

We identified bioregions based on benthic forage community composition using k-means cluster analysis (Legendre et al. 2012) and evaluated change in their spatial extent in the SEBS over time. We focused on the SEBS again due to the availability of diet data over a greater number of years. To perform the analysis, we first standardized composite densities across forage groups by rescaling values to have a maximum of 1 (values were divided by the maximum composite density of each respective forage group). Next, we square-root transformed values to downweight the relative importance of locations and time periods with unusually high densities (Legendre and Legendre 2012). To identify regions and time periods with similar community composition, we performed k-means cluster analysis on the transformed site × year forage community densities. The clustering method is robust to species with variable or patchy distributions and is useful for grouping geographic regions based on species distribution data (Pata et al. 2024). Specifically, we assumed a range of clusters (k = 2 to 15) and selected the cluster number at which additional clusters did not appreciably decrease the total within-cluster sum of squares (Legendre et al. 2012). We subsequently assigned bioregions across all site × year combinations according to cluster classification and examined change in spatial extent over time and differences in community composition between bioregions using principal components analysis (PCA).

Results

Composite measures of forage group density that averaged foraging rates of predators in the model (D_Eq,s,t) and weighted averages based on the relative prevalence of the forage group in predator diets (D_Wt,s,t) were highly correlated with each other and the third density measure based on the reference predator alone (D_R,s,t; Supporting information). Within all forage groups, correlations between density measures was high, exceeding 0.83 with a mean of 0.96 (Supporting information). For simplicity, we therefore limit all subsequent analyses to D_Wt,s,t.

Spatial correlations between survey- and diet-based forage densities were positive on average, but tended to be lower for polychaetes and bivalves (ρ = 0.45 and 0.20, respectively) relative to those for snow and Tanner crabs (0.69 and 0.61, respectively; Fig. 2A). Spatial correlations also varied greatly among years for polychaetes and bivalves, with negative correlations occurring in 12 and 30% of years, respectively (Fig. 2A). In contrast, for the more well-sampled snow and Tanner crabs, spatial correlations were positive in all years and more tightly distributed around the multiyear average value (Fig. 2A). For polychaetes and bivalves, the validation data set (based on direct benthic grabs) was significantly more limited in time and space and thus likely subject to higher sampling error relative to that used for the crab validation (the annual bottom trawl survey). Based on time-averaged maps for snow crab, the diet-based density estimates captured the northerly distribution evident in the survey-based estimates, although the former indicated higher densities over a wider area (Fig. 2B). For Tanner crab, diet-based density estimates were also able to capture spatial patterns evident in the survey-based data, namely, densities concentrated in the southwestern region of the survey area (Fig. 2B).

[IMAGE OMITTED. SEE PDF]

To summarize and highlight differences in the distribution of each forage group, we obtained time-averaged maps of relative densities (Fig. 3). In general, forage groups with higher densities in deeper, offshore waters included brittle stars and Tanner crabs while benthic shrimps were concentrated in shallower, inshore waters (Fig. 3). Benthic amphipods were concentrated mid-shelf, and snow crab also exhibited a strong increase in abundance with latitude (Fig. 3). The remaining forage groups exhibited patterns that were more irregular, for instance, bivalve and miscellaneous crustacean densities were roughly u-shape in distribution, with lower densities mid-shelf and to the north (Fig. 3). These general spatial patterns were also evident across years for each forage group, though spatiotemporal variation was also prevalent to varying degrees (Supporting information). For instance, the offshore dominance of brittle stars and inshore dominance of benthic shrimps largely persisted across years (Supporting information), whereas in high abundance years, densities of polychaetes appear to expand southward and shoreward (Supporting information), benthic amphipods expand southward and across the survey domain (Supporting information), and snow crab expand at the southern end of their distribution (Supporting information).

[IMAGE OMITTED. SEE PDF]

Areas designated as abundance hotspots with moderate and high agreement across predators (> 50 and 80% agreement, respectively) tended to coincide with regions with high densities based on the composite density values (Fig. 3). However, for some forage groups such as polychaetes and miscellaneous worms, the area of high hotspot agreement was small, indicating higher spatial variation among predators in diet patterns for those groups (Fig. 3). In contrast, areas of high hotspot agreement were larger for benthic shrimps, miscellaneous crustaceans, and brittle stars, indicating lower spatial variation among predators. In the case of snow and Tanner crabs, only one predator (the widely distributed Pacific cod) was included in their respective models, thus no disagreement is represented in the hotspot demarcation.

Spatial coverage of diet samples for predators of bivalves, brittle stars, and miscellaneous crustaceans was poor after 2001 (Supporting information). Abundance indices for those forage groups are therefore not shown for those years. Overall, benthic species showed a variety of trends and apparent cycles (Fig. 4), which could be fitted as abundance indices in future models. For example, benthic amphipods and miscellaneous worms both went from relatively high to low abundance from 1990 to 2000. Alternatively, polychaetes and Tanner crabs have cycles from high to low to high abundance from 2010 to 2022, which roughly correspond to cold–warm stanzas.

[IMAGE OMITTED. SEE PDF]

We limited the bioregional analysis to six forage groups with the most complete temporal sampling in the data set: polychaetes, benthic amphipods, benthic shrimp, miscellaneous worms, and snow and Tanner crabs. Five bioregions were identified through cluster analysis with distinct community compositions (Fig. 5). On average, bioregion A occurred in inshore waters and was characterized by relatively high densities of benthic shrimps and low crab densities (Fig. 5). Bioregion B occurred in northern mid-shelf waters and contained higher relative densities of polychaetes and both snow and Tanner crabs (Fig. 5). Bioregion C and D both occurred mid-self but the latter was characterized by higher densities of benthic shrimp, miscellaneous worms, and snow crab while densities for all forage groups in the former were near average (Fig. 5). Finally, bioregion E was characterized by higher densities of Tanner crab and occurred in offshore waters (Fig. 5). To illustrate interannual variation in the spatial configuration of the benthic forage community, we provide bioregion maps for 1996, 2010, and 2016 which correspond to average, cold, and warm years in terms of historical bottom temperatures (Fig. 5D). Overall, bioregions occupying the mid-shelf (B, C, and D) exhibited the largest changes in aerial coverage over time (Fig. 5C). For instance, bioregion D (benthic amphipods, miscellaneous worms, and snow crab) covered ~ 25% of the SEBS from 1984 to 2015, then decreased to ~ 10% thereafter (Fig. 5C). In contrast, bioregion B (polychaetes, snow and Tanner crab) fluctuated in area between 4 and 22% prior to 2015, but showed an increasing trend thereafter, approaching 35% by 2023. Similarly, bioregion C (average relative densities) also fluctuated in the first half of the time series, but exhibited a decreasing trend in the latter half of the time series (Fig. 5C). Long-term trends were less apparent for inshore bioregion A (benthic amphipods) and offshore bioregion E (Tanner crab; Fig. 5C).

[IMAGE OMITTED. SEE PDF]

Discussion

Fish productivity in the EBS strongly depends on benthic invertebrate resources but capacity to monitor infauna and small epibenthos at the ecosystem scale through conventional direct sampling is currently limited. Our analysis helps to address this gap and demonstrates application of a spatiotemporal modeling framework that permits synthesis of diet data from multiple predators that can enable a variety of products including estimation of distribution maps, regional biomass indices, and assessment of forage community patterns. The approach thus offers a potential way to build monitoring capacity for forage groups using diet data that are already routinely collected and could be usefully applied to other systems where similar sampling efforts exist.

Importantly, we show that community structure, and thus foraging opportunities for benthic predators, varies across the shelf and is dynamic. In general, densities of polychaetes, benthic amphipods, and snow crab increase with latitude, brittle star and tanner crab increase offshore, benthic shrimps increase inshore, and spatial patterns for the remaining groups are more irregular. Against these average spatial patterns, we show that interannual variation in abundance also fluctuates widely for each forage group along with distributional shifts. Other spatial distribution data on benthic forage groups, particularly infauna, are highly limited for the EBS. To our knowledge, the only studies in the past 40 years used benthic grab samples collected opportunistically over the period 2006–2021 to examine the distribution of the dominant infaunal groups (Yeung et al. 2010, Yeung and Yang 2018). The results generally show higher polychaete densities on the inner shelf (≤ 50m deep) – both in the NBS and in the SEBS near the Alaska Peninsula. Amphipod densities also trended higher in the north, while bivalves were more widely distributed, with lower-density areas nearshore north of 60°N and on the middle shelf (50–100 m depth) north of the Peninsula (Yeung and Yang 2018). However, this spatially sparse data set covers only a small portion of the EBS and cannot address temporal variability, thus complicating direct comparison with our results. The present analysis has focused on the synthesis of diet data to generate a descriptive understanding of how benthic forage groups change over space and time, but an important next step is to explain how environmental conditions govern variability. We recommend extending the spatiotemporal modeling framework presented here and estimating the effect of candidate environmental covariates or, alternatively, by applying time series analyses (Thorsen et al. 2024b) to abundance indices generated using the model.

The multi-predator modeling approach we employed builds on previous single-predator efforts that adopted a similar spatiotemporal modeling framework (Ng et al. 2021, Grüss et al. 2023), but affords at least three principal benefits. First, including multiple predators with different distributions increases the area sampled by predators in aggregate (and therefore the spatial coverage of stomach samples as well) and is more likely to encompass the forage species or group distribution. As overlap increases between predator and forage species, the resulting diet-based abundance index becomes more informative (Ng et al. 2021). Second, by including stomach content data from multiple predator species in the model, higher frequencies of samples over time and densities over space can be achieved, permitting estimation of finer-resolution spatiotemporal patterns. Last, the precision of parameter estimates and derived quantities may be better in multi-predator relative to single-predator models due to the larger quantity of data used in their fitting (Grüss et al. 2023), while also allowing us to characterize aspects of variation in forage density maps among predators (e.g. hotspot agreement across predators, Fig. 3). That said, inclusion of multiple predators requires careful consideration of how foraging rates are combined to develop composite measures of local forage densities. Other differences in predator behavior, namely, foraging range, may also complicate efforts to develop abundance indices on finer spatial scales – the level of correspondence between stomach content data and local forage densities may differ across predators if foraging ranges also differs. Although the two approaches that we considered were highly correlated, we do not presume this will necessarily be the case in other systems or for other prey types. We encourage additional research into other potential methods for integrating predictions across predators, for example, if sufficient survey-based density estimates are available, data-weighted ensemble approaches or cross-validation optimization methods may be advantageous (Dormann et al. 2018). Further, depending on research objectives, it may also be reasonable to pool diet data, disregarding predator species altogether (Gaichas et al. 2024).

The present analysis demonstrates methods for inferring benthic forage abundances from diet data, but we note that the derived products have several potential applications in the context of research and ecosystem-based fisheries management in the EBS. Foremost, biomass indices for forage species may help to characterize system-level changes in benthic productivity in the EBS and the favorability of conditions for fish growth and potentially survival (Yeung et al. 2010, Yeung and Yang 2014, 2014). Such insights can help provide additional qualitative context for decision-making with respect to total allowable catch limits (Zador et al. 2017, Barbeaux et al. 2020). Regional maps of benthic forage can inform fisheries management decisions regarding spatial management. For example, fisheries managers in the USA are required to identify the location of prey species' habitat in their fisheries management plans (i.e. Component-7 of essential fish habitat, see US Code of Federal Regulations Title-50 § 600.815, ), and this information could be used to minimize fishing impacts on prey abundance. From a modeling perspective, abundance indices could prove useful as covariates in single- and multispecies stock assessments (Ianelli et al. 2016, Holsman et al. 2020, Gaichas et al. 2024) and aid calibration and parameterization of ecosystem models (Reum et al. 2020, Whitehouse et al. 2021), and distribution patterns can help refine spatially explicit EBS models including biogeochemical models coupled to regional ocean models (Gibson and Spitz 2011, Hermann et al. 2016, Kearney et al. 2020). Importantly, the modeling approach employed here, along with recent extensions (e.g. dynamic structural equation models, Thorson et al. 2024b), also provides a framework for evaluating how habitat variables, predation pressure, and disturbances such as those arising from bottom trawling, may drive observed spatiotemporal patterns (Thorson 2019), and are the foci of ongoing research.

Our analysis demonstrated correspondence between diet- and survey-based forage density estimates, which offers encouraging support for inferring abundance patterns using diet data. However, we also note the lower mean and larger spread of annual spatial correlation coefficients for polychaetes and bivalves relative to snow and Tanner crabs. For the former forage groups, weaker and more variable correlations were likely driven by higher sampling error, but modeling artifacts and predator foraging behaviors may also have played important roles. Relative to the crab validation data set, higher levels of sampling error may be driven by: 1) the smaller sampling area of benthic grabs relative to the bottom trawl (0.1 m² versus ~ 0.25 km²); 2) the sampling of fewer stations (23 versus ~ 350 stations each year on average for polychaetes, bivalves, and crabs, respectively); and 3) limited temporal replication of sampling at stations – if spatial correlations differ between subregions, this will impart additional variability. Further, it is possible that the coarse taxonomic levels we aggregated the benthic grab data to (to maximize comparability with the most diet samples) weakened relationships that may exist if finer taxonomic levels are adopted (Yeung et al. 2010).

However, challenges more general to our modeling approach may also have weakened correlations. For instance, the models' average diet information at scales related to knot spacing (here, a minimum distance of 60 km), and it is possible that models fit at finer scales may improve correlations at the cost of computational efficiency. Similarly, predator behavior may limit the potential correlation between diet and local prey densities if, for example, predators forage over large areas or exhibit behaviors such as prey switching in the presence of alternative prey. While we are presently unable to disentangle the relative importance of these different processes, the stronger, more stable correlations associated with the larger validation data set for Tanner and snow crabs indicates that the approach holds merit. That said, we encourage large-scale paired sampling of predator stomachs and macrobenthos to further refine and validate our method for EBS predators.

We envision several lines for continued research regarding spatiotemporal models for stomach content analysis:

• 1) Prey switching: Changes in thinning rates (i.e. attack and capture probability) are confounded with prey or forage densities. These potential changes could be examined by adding covariates to the model that account for the biomass of alternative prey in predator stomachs, which would essentially correct for decreased predator attack probabilities when alternative prey are present. In particular, these models are amenable to joint models including multiple prey and predators, e.g. using ‘tinyVAST'.

• 2) Nonlinear functional responses: Predators may forage over a wider area when prey densities are low, and this would result in a correlation between thinning rates and prey density. Just as thinning rates are analogous to gear catchability in fishery CPUE standardization, this nonlinearity is analogous to hyperstable CPUE (Hilborn and Walters 1992, Wilberg et al. 2009). We imagine that hyperstability would be identifiable given auxiliary information, e.g. when intercalibrating stomach samples relative to conventional sampling methods.

• 3) Process research: Predators typically digest prey faster in warmer waters (Andersen 1999) or may change prey search volume (and thus prey encounter rates) and prey attack rates (Englund et al. 2011, Gliwicz et al. 2018). Analysts could potentially account for variation in ‘thinning rates' attributable to temperature-dependent stomach evacuation, attack, or search rates by including temperature as an offset variable, using laboratory measurements of the various rates. We therefore encourage collaborative research with experimental researchers to measure temperature dependencies (and other habitat-related dependencies) in foraging processes.

• 4) Joint modelling: Finally, we have relied upon diet data derived from the visual identification and sorting of stomach contents from predators sampled with bottom trawl. However, some prey species may be more easily identified (or to a higher taxonomic resolution) using environmental DNA (Traugott et al. 2021) and methods for estimating abundance are emerging (Stoeckle et al. 2021, Fukaya et al. 2022). We therefore recommend efforts to integrate eDNA and visual stomach sampling. Additionally, past research has applied model-based intercalibration to jointly model stomach contents and conventional gears (Grüss et al. 2023). For those species where stomach contents and surveys produce correlated indices (e.g. Tanner and snow crabs), we envision that these two data sets could be jointly modeled, thereby estimating forage habitat utilization across a larger spatial and/or seasonal domain. Ultimately, we encourage development of more synthetic research programs to take advantage of evolving spatiotemporal modeling frameworks and recent advances in how stomach content data are conceptualized (Thorson et al. 2022).

The use of predators-as-samplers has grown rapidly in the last two decades as a practical, complementary method for characterizing the distribution and abundance of forage species (Link 2004, Deroba 2018, Ng et al. 2021, Grüss et al. 2023, Gaichas et al. 2024). Our analysis makes significant inroads with respect to the synthesis of diet data from multiple predator species and provides among the most comprehensive overviews of benthic prey distributions of the EBS to date. Given evidence of rapid environmental changes due to warming (Siddon et al. 2020, Carvalho et al. 2021, Hunt et al. 2022), an important next step is to evaluate potential drivers of change in the benthic community within the context of the larger EBS ecosystem. That said, the biomass indices and distribution patterns obtained under the modeling framework are valuable in their own right, and should be relevant for ecosystem monitoring, heightening both the value of diet data and need to continue sampling into the future.

Acknowledgements

– We would like to thank A. Whitehouse and three anonymous reviewers for helpful comments on early drafts of the manuscript, and past and present members of the Alaska Fisheries Science Center Food Habits Lab for their enthusiastic collection and processing of stomach content data. The findings and conclusions in this article are those of the authors and do not necessarily represent the views of the National Marine Fisheries Service, NOAA. Reference to trade names does not imply endorsement by the National Marine Fisheries Service, NOAA.

Funding

– All authors were supported by NOAA Fisheries.

Author contributions

Jonathan C. P. Reum: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Methodology (equal); Project administration (equal); Visualization (lead); Writing – original draft (lead). James T. Thorson: Formal analysis (equal); Methodology (equal); Writing – original draft (equal). Cynthia Yeung: Data curation (supporting); Writing – review and editing (supporting). Kerim Aydin: Data curation (supporting); Project administration (supporting); Writing – review and editing (supporting).

Data availability statement

Data are available from the figshare digital repository: (Reum et al. 2025).

References

Andersen, N. G. 1999. The effects of predator size, temperature, and prey characteristics on gastric evacuation in whiting. – J. Fish Biol. 54: 287–301.

Andersen, N. G., Chabot, D. and Couturier, C. S. 2016. Modelling gastric evacuation in gadoids feeding on crustaceans. – J. Fish Biol. 88: 1886–1903.

Aydin, K. Y., Gaichas, S., Ortiz, I., Kinzey, D. and Friday, N. 2007. A comparison of the Bering Sea, Gulf of Alaska, and Aleutian Islands large marine ecosystems through food web modeling. – US Dep. Commer., NOAA Tech. Memo. NMFS‐AFSC‐178.

Barbeaux, S. J., Holsman, K. and Zador, S. 2020. Marine heatwave stress test of ecosystem‐based fisheries management in the Gulf of Alaska Pacific cod fishery. – Front. Mar. Sci. 7: 703.

Bilyard, G. R. 1987. The value of benthic infauna in marine pollution monitoring studies. – Mar. Pollut. Bull. 18: 581–585.

Carey, J. M. and Keough, M. J. 2002. Compositing and subsampling to reduce costs and improve power in benthic infaunal monitoring programs. – Estuaries 25: 1053–1061.

Carvalho, K. S., Smith, T. E. and Wang, S. 2021. Bering Sea marine heatwaves: patterns, trends and connections with the Arctic. – J. Hydrol. 600: 126462.

Coyle, K. O., Konar, B., Blanchard, A., Highsmith, R. C., Carroll, J., Carroll, M., Denisenko, S. G. and Sirenko, B. I. 2007. Potential effects of temperature on the benthic infaunal community on the southeastern Bering Sea shelf: possible impacts of climate change. – Deep Sea Res. II 54: 2885–2905.

Deroba, J. J. 2018. Sources of variation in stomach contents of predators of Atlantic herring in the northwest Atlantic during 1973–2014. – ICES J. Mar. Sci. 75: 1439–1450.

Dormann, C. F. et al. 2018. Model averaging in ecology: a review of Bayesian, information‐theoretic, and tactical approaches for predictive inference. – Ecol. Monogr. 88: 485–504.

Edwards, A. M. and Auger‐Méthé, M. 2019. Some guidance on using mathematical notation in ecology. – Methods Ecol. Evol. 10: 92–99.

Englund, G., Öhlund, G., Hein, C. L. and Diehl, S. 2011. Temperature dependence of the functional response. – Ecol. Lett. 14: 914–921.

Fukaya, K., Murakami, H., Yoon, S., Minami, K., Osada, Y., Yamamoto, S., Masude, R., Kasai, A., Miyashita, K., Minamoto, T. and Kondoh, M. 2022. Estimating fish population abundance by integrating quantitative data on environmental DNA and hydrodynamic modelling. – Mol. Ecol. 30: 3067–3067.

Gaichas, S. K., Gartland, J., Smith, B. E., Wood, A. D., Ng, E. L., Celestino, M., Drew, K., Tyrell, A. S. and Thorson, J. T. 2024. Assessing small pelagic fish trends in space and time using piscivore stomach contents. – Can. J. Fish. Aquat. Sci. 81: 990–1012.

Gibson, G. A. and Spitz, Y. H. 2011. Impacts of biological parameterization, initial conditions, and environmental forcing on parameter sensitivity and uncertainty in a marine ecosystem model for the Bering Sea. – J. Mar. Syst. 88: 214–231.

Gliwicz, Z. M., Babkiewicz, E., Kumar, R., Kunjiappan, S. and Leniowski, K. 2018. Warming increases the number of apparent prey in reaction field volume of zooplanktivorous fish. – Limnol. Oceanogr. 63: S30–S43.

Grebmeier, J. M., McRoy, C. P. and Feder, H. M. 1988. Pelagic–benthic coupling on the shelf of the northern Bering and Chukchi seas. 1. Food supply source and benthic biomass. – Mar. Ecol. Prog. Ser. 48: 57–67.

Grebmeier, J. M., Overland, J. E., Moore, S. E., Farley, E. V., Carmack, E. C., Cooper, L. W., Frey, K. E., Helle, J. H., McLaughlin, F. A. and McNutt, S. L. 2006. A major ecosystem shift in the northern Bering Sea. – Science 311: 1461–1464.

Grebmeier, J. M., Bluhm, B. A., Cooper, L. W., Danielson, S. L., Arrigo, K. R., Blanchard, A. L., Clarke, J. T., Day, R. H., Frey, K. E., Gradinger, R. R., Kędra, M., Konar, B., Kuletz, K. J., Lee, S. H., Lovvorn, J. R., Norcross, B. L. and Okkonen, S. R. 2015. Ecosystem characteristics and processes facilitating persistent macrobenthic biomass hotspots and associated benthivory in the Pacific Arctic. – Prog. Oceanogr. 136: 92–114.

Grüss, A., Thorson, J. T., Stawitz, C. C., Reum, J. C. P., Rohan, S. K. and 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). – Prog. Oceanogr. 194: 102569.

Grüss, A., Thorson, J. T., Anderson, O. F., O'Driscoll, R. L., Heller‐Shipley, M. and Goodman, S. 2023. Spatially varying catchability for integrating research survey data with other data sources: case studies involving observer samples, industry‐cooperative surveys, and predators as samplers. – Can. J. Fish. Aquat. Sci. 80: 1595–1615.

Haflinger, K. 1981. A survey of benthic infaunal communities of the southeastern Bering Sea shelf. – In: Hood, D. W. and Calder, J. A. (eds), The eastern Bering Sea shelf: oceanography and resources. US Department of Commerce, NOAA, pp. 1091–1104.

Hale, S. S., Buffum, H. W., Kiddon, J. A. and Hughes, M. M. 2017. Subtidal benthic invertebrates shifting northward along the US Atlantic Coast. – Estuaries Coast. 40: 1744–1756.

Hartig, F. and Lohse, L. 2020. DHARMa: residual diagnostics for hierarchical (multi‐level/mixed) regression models, ver. 0.3. 3.0. – https://cran.r‐project.org/web/packages/DHARMa/vignettes/DHARMa.html.

Hermann, A. J., Gibson, G. A., Bond, N. A., Curchitser, E. N., Hedstrom, K., Cheng, W., Wang, M., Cokelet, E. D., Stabeno, P. J. and Aydin, K. 2016. Projected future biophysical states of the Bering Sea. – Deep Sea Res. II 134: 30–47.

Hilborn, R. and Walters, C. J. 1992. Quantitative fisheries stock assessment: choice, dynamics, and uncertainty. – Chapman and Hall.

Holsman, K. K., Haynie, A. C., Hollowed, A. B., Reum, J. C. P., Aydin, K., Hermann, A. J., Cheng, W., Faig, A., Ianelli, J. N. and Kearney, K. A. 2020. Ecosystem‐based fisheries management forestalls climate‐driven collapse. – Nat. Commun. 11: 4579.

Hunt, G. L., Yasumiishi, E. M., Eisner, L. B., Stabeno, P. J. and Decker, M. B. 2022. Climate warming and the loss of sea ice: the impact of sea‐ice variability on the southeastern Bering Sea pelagic ecosystem. – ICES J. Mar. Sci. 79: 937–953.

Ianelli, J., Holsman, K. K., Punt, A. E. and Aydin, K. 2016. Multi‐model inference for incorporating trophic and climate uncertainty into stock assessments. – Deep Sea Res. II 134: 379–389.

Jewett, S. C. and Feder, H. M. 1981. Epifaunal invertebrates of the continental shelf of the eastern Bering and Chukchi Seas. – East Bering Sea Shelf Oceanogr. Resour. 2: 1131–1153.

Jørgensen, L. L., Renaud, P. E. and Cochrane, S. K. J. 2011. Improving benthic monitoring by combining trawl and grab surveys. – Mar. Pollut. Bull. 62: 1183–1190.

Kass, R. E. and Steffey, D. 1989. Approximate Bayesian inference in conditionally independent hierarchical models (parametric empirical Bayes models). – J. Am. Stat. Assoc. 84: 717–726.

Kearney, K., Hermann, A., Cheng, W., Ortiz, I. and Aydin, K. 2020. A coupled pelagic–benthic–sympagic biogeochemical model for the Bering Sea: documentation and validation of the BESTNPZ model (v2019. 08.23) within a high‐resolution regional ocean model. – Geosci. Model Dev. 13: 597–650.

Kotwicki, S., Lauth, R. R., Williams, K. and Goodman, S. E. 2017. Selectivity ratio: a useful tool for comparing size selectivity of multiple survey gears. – Fish. Res. 191: 76–86.

Kristensen, K., Nielsen, A., Berg, C. W., Skaug, H. and Bell, B. 2016. TMB: automatic differentiation and Laplace approximation. – J. Stat. Softw. 70: 1–21.

Lang, G. M., Derrah, C. W. and Livingston, P. A. 2003. Groundfish food habits and predation on commercially important prey species in the eastern Bering Sea from 1993 through 1996. – NOAA, National Marine Fisheries Service, Alaska Fisheries Science Center.

Legendre, P. and Legendre, L. 2012. Numerical ecology. – Elsevier.

Lindgren, F., Rue, H. and Lindström, J. 2011. An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. – J. R. Stat. Soc. B 73: 423–498.

Link, J. S. 2004. Using fish stomachs as samplers of the benthos: integrating long‐term and broad scales. – Mar. Ecol. Prog. Ser. 269: 265–275.

Livingston, P. A., Aydin, K., Buckley, T. W., Lang, G. M., Yang, M.‐S. and Miller, B. S. 2017. Quantifying food web interactions in the North Pacific – a data‐based approach. – Environ. Biol. Fishes 100: 443–470.

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

McDonald, J., Feder, M. and Hoberg, M. 1981. Bivalve mollusks of the southeastern Bering Sea. – In: Hood, D. W. and Calder, J. A. (eds), The eastern Bering Sea shelf oceanography and resources. US Department of Commerce, NOAA.

Moore, S. E. and Grebmeier, J. M. 2018. The distributed biological observatory: linking physics to biology in the Pacific Arctic region. – Arctic 71(Suppl. 1) 1: 1–7.

Ng, E. L., Deroba, J. J., Essington, T. E., Grüss, A., Smith, B. E. and Thorson, J. T. 2021. Predator stomach contents can provide accurate indices of prey biomass. – ICES J. Mar. Sci. 78: 1146–1159.

Pata, P. R., Galbraith, M., Young, K., Margolin, A. R., Perry, I. R. and Hunt, B. P. V. 2024. Data‐driven determination of zooplankton bioregions and robustness analysis. – MethodsX 12: 1146.

Reum, J. C. P., Blanchard, J. L., Holsman, K. K. and Aydin, K. 2019. Species‐specific ontogenetic diet shifts attenuate trophic cascades and lengthen food chains in an exploited ecosystem. – Oikos 128: 1051–1064.

Reum, J. C. P., Blanchard, J. L., Holsman, K. K., Aydin, K., Hollowed, A. B., Hermann, A. J., Cheng, W., Faig, A., Haynie, A. C. and Punt, A. E. 2020. Ensemble projections of future climate change impacts on the eastern Bering Sea food web using a multispecies size spectrum model. – Front. Mar. Sci. 7: 124.

Reum, J. C. P., Thorson, J., Yeung, C. and Aydin, K. 2025. Data from: Assessing benthos through predator stomach contents: spatiotemporal modeling of abundance and habitat use. – figshare digital repository, 10.6084/m9.figshare.28941317.

Reynoldson, T. B. and Metcalfe‐Smith, J. L. 1992. An overview of the assessment of aquatic ecosystem health using benthic invertebrates. – J. Aquat. Ecosyst. Health 1: 295–308.

Siddon, E. 2024. Ecosystem status report 2024: eastern Bering Sea, stock assessment and fishery evaluation report. – North Pacific Fishery Management Council.

Siddon, E. C., Zador, S. G. and Hunt Jr, G. L. 2020. Ecological responses to climate perturbations and minimal sea ice in the northern Bering Sea. – Deep Sea Res. II 181–182: 104914.

Smith, J. Q. 1985. Diagnostic checks of non‐standard time series models. – J. Forecast 4: 283–291.

Stevenson, D. E. and Lauth, R. R. 2012. Latitudinal trends and temporal shifts in the catch composition of bottom trawls conducted on the eastern Bering Sea shelf. – Deep Sea Res. II 65–70: 251–259.

Stoeckle, M. Y., Adolf, J., Charlop‐Powers, Z., Dunton, K. J., Hinks, G. and VanMorter, S. M. 2021. Trawl and eDNA assessment of marine fish diversity, seasonality, and relative abundance in coastal New Jersey, USA. – ICES J. Mar. Sci. 78: 293–304.

Stoker, S. W. 1978. Benthic invertebrate macrofauna of the Eastern continental shelf of the Bering and Chukchi Seas. – PhD thesis, Univ. of Alaska, USA.

Stoker, S. W. 1981. Benthic invertebrate macrofauna of the eastern Bering. – In: Hood, D. W. and Calder, J. A. (eds), The eastern Bering Sea shelf: oceanography and resources. US Department of Commerce, NOAA, pp. 1069–1090.

Thorson, J. T. 2019. Guidance for decisions using the vector autoregressive spatio‐temporal (VAST) package in stock, ecosystem, habitat and climate assessments. – Fish. Res. 210: 143–161.

Thorson, J. T. and Kristensen, K. 2016. Implementing a generic method for bias correction in statistical models using random effects, with spatial and population dynamics examples. – Fish. Res. 175: 66–74.

Thorson, J. T., Shelton, A. O., Ward, E. J. and Skaug, H. J. 2015. Geostatistical delta‐generalized linear mixed models improve precision for estimated abundance indices for west coast groundfishes. – ICES J. Mar. Sci. 72: 1297–1310.

Thorson, J. T., Arimitsu, M. L., Levi, T. and Roffler, G. H. 2022. Diet analysis using generalized linear models derived from foraging processes using R package mvtweedie. – Ecology 103: e3637.

Thorson, J. T., Barnes, C. L., Friedman, S. T., Morano, J. L. and Siple, M. C. 2023. Spatially varying coefficients can improve parsimony and descriptive power for species distribution models. – Ecography 2023: e06510.

Thorson, J. T., Anderson, S. C., Goddard, P. and Rooper, C. N. 2024a. tinyVAST: R package with an expressive interface to specify lagged and simultaneous effects in multivariate spatio‐temporal models. – Global Ecol. Biogeogr. 34: e70035.

Thorson, J. T., Andrews, A. G., Essington, T. E. and Large, S. I. 2024b. Dynamic structural equation models synthesize ecosystem dynamics constrained by ecological mechanisms. – Methods Ecol. Evol. 15: 744–755.

Traugott, M., Thalinger, B., Wallinger, C. and Sint, D. 2021. Fish as predators and prey: DNA‐based assessment of their role in food webs. – J. Fish Biol. 98: 367–382.

Warton, D. I., Thibaut, L. and Wang, Y. A. 2017. The PIT‐trap – a “model‐free” bootstrap procedure for inference about regression models with discrete, multivariate responses. – PLoS One 12: e0181790.

Whitehouse, G. A., Aydin, K. Y., Hollowed, A. B., Holsman, K. K., Cheng, W., Faig, A., Haynie, A. C., Hermann, A. J., Kearney, K. A., Punt, A. E. and Essington, T. E. 2021. Bottom–up impacts of forecasted climate change on the eastern Bering Sea food web. – Front. Mar. Sci. 8: 624301.

Wilberg, M. J., Thorson, J. T., Linton, B. C. and Berkson, J. 2009. Incorporating time‐varying catchability into population dynamic stock assessment models. – Rev. Fish. Sci. 18: 7–24.

Yeung, C. and McConnaughey, R. A. 2006. Community structure of eastern Bering Sea epibenthic invertebrates from summer bottom‐trawl surveys 1982 to 2002. – Mar. Ecol. Prog. Ser. 318: 47–62.

Yeung, C. and Yang, M.‐S. 2014. Habitat and infauna prey availability for flatfishes in the northern Bering Sea. – Polar Biol. 37: 1769–1784.

Yeung, C. and Yang, M.‐S. 2017. Habitat quality of the coastal southeastern Bering Sea for juvenile flatfishes from the relationships between diet, body condition and prey availability. – J. Sea Res. 119: 17–27.

Yeung, C. and Yang, M.‐S. 2018. Spatial variation in habitat quality for juvenile flatfish in the southeastern Bering Sea and its implications for productivity in a warming ecosystem. – J. Sea Res. 139: 62–72.

Yeung, C. and Yang, M. S. 2025.Macroinfauna communities on the Bering Sea continental shelf. – Mar. Ecol. Prog. Ser. 759: https://doi.org/10.3354/meps14841.

Yeung, C., Yang, M.‐S. and McConnaughey, R. A. 2010. Polychaete assemblages in the south‐eastern Bering Sea: linkage with groundfish distribution and diet. – J. Mar. Biol. Assoc. UK 90: 903–917.

Zador, S. G., Holsman, K. K., Aydin, K. Y. and Gaichas, S. K. 2017. Ecosystem considerations in Alaska: the value of qualitative assessments. – ICES J. Mar. Sci. 74: 421–430.