Weak effects of local prey density and spatial overlap on predation intensity in a temperate marine ecosystem
Max Lindmark, Christopher A. Griffiths, Valerio Bartolino, Viktor Thunell, Federico Maioli, Sean C. Anderson, Andrea Belgrano, Michele Casini, Katarzyna Nadolna‐Ałtyn, Joanna Pawlak, Marzenna Pachur, Marcin Rakowski, Karolina Wikström, Murray S. A. Thompson, Mayya Gogina
TL;DR
This study finds weak links between prey availability and predation by cod in the Baltic Sea, challenging assumptions in marine ecosystem models.
Contribution
The study uses spatiotemporal models and long-term data to assess predator-prey interactions in a temperate marine ecosystem.
Findings
Predator-prey overlap has decreased over three decades in the Baltic Sea.
Only Saduria shows a clear link between prey availability and predation by cod.
Pelagic prey like herring and sprat show weak connections to cod predation patterns.
Abstract
Quantifying the impact of lower trophic level species abundance on higher trophic level predators (and vice versa) is critical for understanding marine ecosystem dynamics and for implementing ecosystem‐based management. Trophic ecosystem models generally predict a tight coupling between prey and fish predators, such that higher abundance of lower trophic species increases the abundance of higher trophic level predators. This assumes that predator feeding rates are limited by prey availability to some degree. Despite being a key component of predator–prey interactions and multispecies fisheries management, relatively few studies have assessed the impacts of prey availability on predation patterns of mobile, generalist fish predators using spatiotemporal models and local‐scale stomach content, predator, and prey data. In this study, we explore the association between local density of key…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
FIGURE 1
FIGURE 2
FIGURE 3
FIGURE 4
FIGURE 5
FIGURE 6
FIGURE 7
FIGURE 8
FIGURE 9| Prey | ΔAIC | ||
|---|---|---|---|
| Breakpoint | Linear | No prey covariate | |
| Herring | 0 | 274 | 0.09 |
|
| 0 | 38.7 | 51.5 |
| Sprat | 1.78 | 2616 | 0 |
- —German Federal Ministry of Education and Research
- —Oscar and Lili Lamm Memorial Foundation10.13039/100007794
- —Natural Environment Research Council10.13039/501100000270
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsMarine and fisheries research · Isotope Analysis in Ecology · Marine and coastal ecosystems
INTRODUCTION
Temperate marine ecosystems are shaped by complex combinations of bottom‐up and top‐down processes (Lynam et al., 2017). In coastal or upwelling systems, these two processes can be linked at intermediate trophic levels by small pelagic planktivorous fish, which can control lower trophic levels with top‐down effects and upper trophic levels via bottom‐up effects—also referred to as wasp‐waist control (Bakun, 2006; Cury et al., 2000). Ecological theory generally suggests a tight coupling between prey and predators (Poggiale, 1998) such that higher prey densities benefit predator populations. This positive link between prey abundance and predator performance has been supported by a range of trophic ecosystem models (Chagaris et al., 2020; Pikitch et al., 2014; Smith et al., 2011), although this may vary depending on structural model assumptions (Walters et al., 2016). Such findings have led in some circumstances to management strategies that aim to reduce fishing pressure on lower trophic level fish to avoid negative impacts on higher trophic levels (Chagaris et al., 2020; Cury et al., 2011; Pikitch et al., 2012).
Empirical support for this form of bottom‐up or “donor effect” (i.e., positive effects of prey on the predator but no clear negative effects of predators on prey, sensu Jennings & Kaiser, 1998) from intermediate to upper trophic levels is, however, fairly limited (Jensen et al., 2012). This could be because marine predatory fish are typically mobile generalists, with diets changing through their ontogeny, which implies many but weak trophic links (Jennings & Kaiser, 1998; Strong, 1992). Not surprisingly, one of the clearer examples of a donor effect is the decline in productivity of Atlantic cod (Gadus morhua) in the Barents Sea. This decline occurred because cod had limited ability to switch prey due to low diversity in the intermediate trophic levels, leading to reduced productivity after the main prey species, capelin (Mallotus villosus) and herring (Clupea harengus), collapsed (Hamre, 1994; Jensen et al., 2012; Mehl & Sunnanå, 1991).
Another reason why these effects are difficult to detect relates to spatial dynamics and the scale at which trophic interactions occur (Hunsicker et al., 2011). Often trophic interactions are assessed by relating time series of population abundance of predators and prey (Hilborn et al., 2017; Jennings & Kaiser, 1998; Overholtz & Link, 2007), which neglects the important spatial dimension. While there are several examples of the importance of local scale prey availability and donor effects on upper trophic levels, most are for central place foragers (e.g., for seabirds, Crawford, 2007; Cury et al., 2011; Hentati‐Sundberg et al., 2021; Robinson et al., 2015), and few involve commercially exploited mobile fish predators (Hilborn et al., 2017; Pikitch et al., 2018) (but see Fall et al., 2021). Recent studies have illustrated the potential of spatiotemporal modeling of predator and prey density and stomach content data to quantify fine‐scale spatial and temporal variability in diet and to derive model‐based indices of consumption (Gartland & Latour, 2024; Goodman et al., 2022; Grüss et al., 2020). For example, Gartland and Latour (2024) related model‐based consumption indices derived from spatiotemporal models to stock‐level biomass of prey and found positive associations. In a similar study, Goodman et al. (2022) correlated annual spatial overlap and predation indices, and found that the support for such correlations varied across species. These studies did not, however, include spatially explicit prey covariates, which may be important, because local prey densities can be high and sufficient for predators even though total prey population abundance is low (Hilborn et al., 2017), and because the spatial overlap between predators and their prey can vary in the studied time period.
In this study, we investigated the relationship between local prey availability and the relative mass of prey in predator stomachs, as well as the relationship between spatial overlap and predation at a population level, using the Baltic Sea as a case study. We used three decades of biomass and stomach content data from surveys for the predator Atlantic cod (G. morhua), and biomass data for some of its main prey, sprat (Sprattus sprattus), herring (Clupea harengus), and the benthic isopod Saduria entomon (henceforth only Saduria). The central Baltic Sea ecosystem is a species‐poor ecosystem where sprat and herring make up more than 75% of the diet by mass in cod around 35 cm (Kulatska et al., 2019; Lindmark et al., 2025). Moreover, within this time period, the feeding rates on sprat and Saduria, as well as the growth, condition, and size‐at‐maturity of cod have declined substantially (Lindmark et al., 2023; Mion et al., 2021; Neuenfeldt et al., 2020; Svedäng et al., 2024). Together with an increased natural mortality (Eero et al., 2023; International Council for the Exploration of the Sea [ICES], 2022a), this has severely impacted the conservation status of cod, to the degree that even in the absence of a commercial targeted fishery (by European Union [EU] countries), the stock is expected to remain below its biological limit reference point in the near future (ICES, 2021). Hence, it is critically important to understand the drivers behind these changes in the physiological performance of cod in the Baltic Sea. Reduced feeding opportunities on benthic and pelagic prey have been proposed as two of the possible underlying causes (Casini et al., 2016; Neuenfeldt et al., 2020). Reduced feeding on benthic prey is thought to be due to increased competition for dwindling benthic prey resources mainly linked to the deterioration of the benthic habitats (Neuenfeldt et al., 2020), while reduced feeding on pelagic prey has been hypothesized to be due to a reduction of the prey in the main areas of cod occurrence (Casini et al., 2016; ICES, 2023). This has led to suggestions that pelagic fisheries should be limited in the current main distribution area of cod to improve growth (Casini et al., 2016; Eero et al., 2012; ICES, 2015). However, it is largely unknown how the overall local‐scale overlap and encounter rates have changed over time and how the local prey availability affects feeding opportunities for cod.
Here, we aim to answer the following research questions: (1) Is the relative mass of sprat, herring, and Saduria in cod stomachs related to the local availability of these prey? (2) Have the relative masses of, and predation on, these prey in the diet changed over space and time? (3) Has the spatial overlap between cod and its prey changed over time? (4) Do spatial overlap indices correlate with population‐level predation indices over time?
METHODS
Data
Stomach data
Stomach content data for Baltic cod have been collected by national institutes (opportunistically or within designated programs) over the past decades. In this study, we use mainly data collated in EU‐funded projects (Huwer et al., 2014; Jacobsen et al., 2023), available at the newly updated ICES stomach content database (https://www.ices.dk/data/data-portals/Pages/Stomach-content.aspx). We complement these data with data collated in Lindmark et al. (2025), and historical data (Huwer et al., 2014). We used 31 years (1993–2023) of data, mainly collected on the biannual Baltic International Trawl Survey (BITS) conducted in the first and fourth quarters of the year, but also from other cruises conducted at other times of the year (see Huwer et al. (2014) and Jacobsen et al. (2023) for a detailed description of the data sources). Therefore, the spatial coverage of stomach sampling has varied over time (see Figure 1 for the spatiotemporal distribution of data).
Location of stomach samples in the southern Baltic Sea. Colored contours correspond to depth, the fill color of the points corresponds to year, and the size (area) of the points corresponds to sample size by year at that location. Rectangles correspond to International Council for the Exploration of the Sea (ICES) rectangles.
Over time, the treatment of regurgitated stomachs (and thereby the classification of empty stomachs) has changed (Neuenfeldt et al., 2020). In the early part of the time series, gall bladder status was used to separate nonfeeding predators from feeding predators with stomachs regurgitated in the trawl. In recent years, recording gall bladder status is not mandatory, and predators with visible signs of regurgitation are replaced on board. These differences in sampling over time can lead to biased estimates of the proportion of empty stomachs (Neuenfeldt et al., 2020). We acknowledge there may be trends due to sampling in the proportion of nonfeeding cod and that we cannot assess this bias. Therefore, we opted to include all stomachs in our main analysis, and reran the analysis with empty stomachs omitted (across all years, the mean proportion of empty stomachs was 0.19).
For each cod predator, we calculated the total mass of sprat, herring and Saduria in the stomachs. Any of these prey were found in 43% of predator stomachs. In cases where individual prey length but not prey mass was available, we estimated mass using species‐specific condition factors and mass‐length exponents. The condition factor and exponent for Saduria are means for isopods in Robinson et al. (2010), and for fish prey (herring and sprat) they were retrieved from FishBase (Froese & Pauly, 2025). These individual‐level prey masses were used to calculate prey‐mass‐per‐predator‐mass—hereafter referred to as “relative prey mass.” The size distribution of cod, and the distribution of relative prey masses can be seen in Appendix S1: Figures S1 and S2. Similarly, cod mass was estimated from cod length if missing in the data. Because the length–mass relationship in Baltic cod has varied substantially over the time period, affecting the relative prey mass, we estimated annual values of the condition factor from the trawl survey BITS (see Biomass density models). After data processing, 33,243 cod stomachs were available for analysis.
Prey data
ICES rectangle‐level (Figure 1) biomass estimates for sprat and herring were acquired from the ICES Baltic International Fish Survey Working Group (WGBIFS) database for the Baltic International Acoustic Survey (BIAS) (https://www.ices.dk/community/groups/pages/WGBIFS.aspx). As in Lindmark et al. (2023), biomass density of Saduria was extracted from a habitat distribution model coupled with modeled hydrographical data from the regional coupled ocean biogeochemical model Ecological Regional Ocean Model (ERGOM) (Gogina et al., 2020; Neumann et al., 2021). The model was trained to the time period 1981–2019 and predicted for the time period 1993–2019 (but note that this prediction is constant over time and therefore represents the core Saduria habitat).
Predator biomass density
To calculate spatially explicit, weighted predation metrics and predator–prey overlap metrics (Figure 2), we modeled the spatiotemporal distribution of cod using catch per unit effort data (CPUE, in numbers per hour) by length class from the fishery‐independent Baltic International Trawl Survey (BITS) conducted in the first and fourth quarter between 1993 and 2023 in the ICES subdivisions 24–28. We used data from the ICES trawl survey database DATRAS (https://www.ices.dk/data/data-portals/Pages/DATRAS.aspx). CPUE data were standardized based on gear dimensions and towing speed (TVL trawl with 75‐m sweeps at three knots) to units of kilograms per square kilometer, following Lindmark et al. (2023) and Orio et al. (2017), using length–mass relationships fitted by year to convert from numbers‐at‐length to mass.
Chart describing the data and modeling process to acquire predation indices (population level and per capita) as well as predator–prey overlap. For the definition of Pp,i and Pc,i, see Equations (10) and (11). Research questions are indicated in darker blue boxes.
Environmental data
We included dissolved sea bottom oxygen concentration (in PSU), sea bottom temperature (in degrees Celsius), and sea bottom salinity (in per mille) from the Copernicus Marine Service Baltic Monitoring and Forecasting Centre (BAL MFC) as covariates in our biomass density models. Dissolved oxygen values at the sea floor stem from the Baltic Sea biogeochemical model (https://doi.org/10.48670/moi-00009), which is based on ERGOM (https://ergom.net/) (Neumann et al., 2021), coupled with the NEMO ocean model (Madec et al., 2023). Sea floor temperature and salinity stem from the Baltic Sea Physical Reanalysis (https://doi.org/10.48670/moi-00013), based on simulations from the NEMO 3D ocean‐ice model version 4.0 (Gurvan et al., 2019). These variables were matched to the survey catch data on a monthly resolution. We also included depth (in meters) in the model, which was extracted from the EMODnet Bathymetry project (https://emodnet.ec.europa.eu/en/bathymetry) (EMODnet Bathymetry Consortium, 2022).
Spatiotemporal modeling framework
Model description
We used spatiotemporal generalized linear mixed‐effects models (GLMMs) to model stomach contents and biomass density of cod. Per‐capita and population‐level predation, and predator–prey overlap, were calculated from predictions onto a 3 × 3 km spatiotemporal grid. Figure 2 illustrates the workflow and how the models and data come together. The full model can be written as follows, but note that the biomass density and diet models do not contain all these terms:
where ys,t is the response variable (relative prey mass or cod biomass density) in location s at time t, μ is the mean, f−1 is the inverse link function, Xmain, Xtvc, and Xsvc are design matrices for fixed‐effects, time‐varying coefficients, and spatially varying coefficients, respectively. Their corresponding coefficient vectors are β, γg, and ζk, where each γg represents a temporally varying coefficient for covariate g, and each ζk represents a spatial field for the k‐th spatially varying coefficient. All fixed‐effect covariates were standardized by subtracting their mean and dividing by their SD. The parameter αg represents a random intercept for month m. The spatially varying coefficients (ζk), spatial (ωs), and spatiotemporal random effects (ϵs,t), are assumed drawn from Gaussian Markov random fields (GMRFs) with covariance matrices (i.e., inverse precision matrices) Σζ, Σω, and Σϵ constrained by Matérn covariance functions (Lindgren et al., 2011; Rue et al., 2009). The spatially varying coefficients ζk are included in the biomass density model to allow the distribution of cod to vary between quarters in addition to the variation given by changes in dynamic covariates, while ωs and ϵs,t reflect spatially correlated latent effects that are constant through time and that vary through time, respectively. The Stochastic Partial Differential Equation (SPDE) approach (Lindgren et al., 2011), which links continuous Gaussian random fields with discretely indexed GMRFs, requires piece‐wise linear basis functions defined by a triangulated mesh. We defined this mesh using triangles with a cutoff distance (minimum distance between vertices) of 8 km and 10 km for diet and biomass models, respectively, and kept all other arguments in fm_rcdt_2d_inla() (fmesher R package; Lindgren, 2023) at their defaults (Appendix S1: Figures S3 and S4). Across all models, the cutoff distance was between 2.5 and 10 times smaller than the estimated Matérn range, which is defined as the distance where spatial correlation effectively disappears (≈0.13, Lindgren et al., 2011).
To propagate uncertainty in both stomach content predictions and cod density predictions when calculating predation and overlap indices (see Predation indices and Spatiotemporal overlap indices), we simulated 500 draws from the joint parameter precision matrix to make model predictions on the grid. For each draw, we calculated overlap or predation metrics and we present the median, mean, or CV of these draws. We predict from the model for a cod length of 33 cm, which is the mean in the diet data.
Diet models
The relative prey mass models were fit to each prey species separately. Since these response variables are positive continuous and contain zeroes, we modeled them with either a Tweedie distribution (herring and sprat) or as a “Poisson‐link” delta–gamma model (Thorson, 2018) (Saduria), which has the flexibility of a classic delta or hurdle model (Aitchison, 1955) but enables a simpler interpretation of covariates due to log links on both linear predictors (Anderson et al., 2024; Thorson, 2018). The model family (delta‐gamma or Tweedie) was chosen based on convergence diagnostics. We included a linear depth effect on all prey, as it has been shown previously to affect cod diets (Pachur & Horbowy, 2013). Month was included as a random effect in the Saduria model, but not the herring and sprat models because the SD of the month random effect was estimated near zero. We included a linear effect of predator length, reflecting the underlying ontogenetic diet shift cod undergo. To be able to interpolate across a missing year (2011), we modeled the intercept as a time‐varying coefficient following an AR(1) process:
where ρ is the correlation between subsequent years and σ2 is the variance. We included spatially correlated latent effects (ωs) that are constant over time and that vary independently each year (ϵs). Lastly, we included prey biomass or biomass densities as covariates. For the sprat and herring models, this was the estimated biomass per ICES rectangle for ages 1–8+, while for Saduria, we extracted rasterized local biomass densities. For each prey, we explored three models: a breakpoint (hockey stick) function (corresponding to a type I functional response with saturation or a type II functional response), a linear effect, and one without prey availability as covariate. In the breakpoint term, prey×b0 is the effect below the threshold, and b0×b1 is the value at the asymptote. We compared the models using the (marginal) Akaike information criterion (AIC; Akaike, 1973).
Biomass density models
We modeled cod biomass density as a Poisson‐link delta–gamma model (Thorson, 2018) and included year and quarter as factors, where the latter was in addition modeled as a spatially varying effect, linear effects of salinity, temperature and temperature squared, as well as a breakpoint effect of oxygen, to reflect that dissolved oxygen tends to correlate with biomass density up to a certain point (Essington et al., 2022). We also included a time‐varying effect of depth and depth squared following a random walk. This was to allow the unimodal depth preferences to change over time (English et al., 2022), in line with the shallowing that has been observed for eastern Baltic cod (Lindmark et al., 2023; Orio et al., 2019).
Lastly we included both spatial random effects that are constant in time (ωs) and spatial random effects that are independent each year (ϵs,t).
Model fitting
We fit the spatiotemporal models with the R version 4.3.2 (R Core Team, 2024) package sdmTMB (Anderson et al., 2024), version 0.6.0.9023. sdmTMB uses automatic differentiation and the Laplace approximation from the R package TMB (Kristensen et al., 2016) and sparse matrix structures to set up the SPDE‐based GMRFs from the R package fmesher (Lindgren, 2023). We estimated parameters via maximum marginal likelihood using the nonlinear minimizer nlminb (R Core Team, 2024). We confirmed that the optimization was consistent with convergence by checking that the Hessian matrix was positive definite and the maximum absolute log‐likelihood gradient with respect to fixed effects was <0.001. We evaluated consistency of the models with the data by calculating simulation‐based randomized quantile residuals (Dunn & Smyth, 1996; Hartig, 2022) (Appendix S1: Figures S5 and S6). When calculating expected values for the purposes of residual calculation, we took a single draw of the random effects from their multivariate normal distribution (Waagepetersen, 2006) rather than using the empirical Bayes random effects estimates. This acknowledges that the random effects are estimated from a distribution (Thorson & Kristensen, 2024, p. 41; Waagepetersen, 2006).
Predation indices
We followed the approach presented in Goodman et al. (2022) to calculate density‐weighted per‐capita and population‐level predation intensity based on model‐predicted relative prey masses and predator densities across the spatiotemporal grid. Population‐level predation intensity, Pp,i, for year i was calculated as
where Ri,j is the prey‐specific relative prey masses (in kilograms per kilogram) in grid cell j, Di,j is the predicted cod density (in kilograms per square kilometer), and Aj is the area of the grid cell (in square kilometers). This represents a spatially explicit, density‐weighted measure of predation intensity (an instantaneous “snapshot” of total mass of a prey species in cod stomachs in units of kilograms) (Goodman et al., 2022). Temporal trends in predation intensity in the domain were acquired by summing grid‐level predictions by year. Prior to calculating grid‐cell‐level predation, we omitted grid cells with cod biomass density predictions greater than the 99.99th percentile across all simulations, and relative prey mass predictions >1 (less than 0.001% of rows in simulated values of relative sprat mass). These filters did not have a qualitative impact on the calculated metrics, but made sure draws where ecologically realistic. We also omitted areas deeper than 130 m in the prediction grid, since trawl surveys are not conducted at those depths.
Since both cod and pelagic species have undergone shifts in their spatial distribution in this time window (Bartolino et al., 2017; ICES, 2023; Lindmark et al., 2023; Orio et al., 2019), and the average feeding rate of cod has declined over time (Neuenfeldt et al., 2020), we also wanted to disentangle the effect of distribution shifts from changes in mean feeding rates. This was done by dividing the population‐level predation (Pp,i) by the population‐level cod biomass in year i (Goodman et al., 2022), which yields a weighted average prey biomass per unit predator biomass (per‐capita predation), Pc,i:
Spatiotemporal overlap indices
Predator–prey overlap metrics were calculated from grid‐level (3 × 3 km) predictions of cod biomass density, biomass density of Saduria, and biomass of sprat and herring at the ICES rectangle level (Figure 1). There are numerous ways to calculate overlap metrics with slightly different interpretations (for a review, see Carroll et al., 2019). In this study, we use the “local index of collocation” overlap metric (Pianka, 1973), as in Goodman et al. (2022):
where overlap in year i is calculated from the proportions of total biomass of cod (predi,j) and its prey (preyi,j) in grid cell j. We use the same grid as for the stomach content predictions. This metric ranges between 0 and 1 and estimates co‐occurrence using correlations between predator and prey densities at the grid scale and is suitable for estimating encounter rates (Carroll et al., 2019; Goodman et al., 2022). When visualizing the overlap in space, we omit the summation across grid cells.
For pelagic prey (sprat and herring), we used only cod predictions in the fourth quarter, since that is when the hydroacoustic survey (BIAS) takes place. After confirming that the difference between quarters was minimal for overlap with Saduria with respect to trends, we presented only results for quarter 4 also for Saduria. Given that the spatial predictions of Saduria densities are constant over time, changes in overlap with Saduria are only driven by changes in the distribution of cod, while in reality, Saduria also likely has changed their spatial distribution (Gogina et al., 2020).
To quantify how overlap between cod and its prey was related to the predation intensity, we tested the correlation between annual predation intensity (per capita and population level) and the annual spatial overlap.
Lastly, to highlight trends over time in per‐capita and population‐level predation and spatial overlap, we fit generalized additive models (GAMs) to the annual indices, with prey‐specific smooth effects of year. We assumed Gamma‐distributed errors and a log link for the predation indices, and Beta‐distributed errors and a logit link for the overlap indices. The models were fitted with the R package brms (Bürkner, 2017). We used default priors, that is, Student‐t3,0,2.5 for the intercept and the SDs of spline coefficients, and flat priors for the prey coefficients. The shape parameter in the Gamma model, and the precision parameter in the Beta model, were given a Gamma(0.01, 0.01) prior. We visualized predictions by summarizing draws from the expectation of the posterior predictive distribution, using the R package tidybayes (Kay, 2023).
RESULTS
The relative mass of a prey species in cod stomachs was positively related to the local density of that prey for Saduria, but not for sprat or herring. For sprat, the model without prey biomass was favored by AIC (Table 1). For herring, the breakpoint model and the model without herring were indistinguishable in terms of AIC (Table 1), and the herring covariate was negatively related to herring in the stomachs of cod. Therefore, further analysis was based on the simpler model without the breakpoint herring covariate. The breakpoint model had the lowest AIC among the Saduria models; however, the conditional predictions showed high uncertainty on the total prediction (both model components combined) (Appendix S1: Figure S7; Table S1).
We found clear spatial patterns in both stomach contents and predation indices (Figures 3 and 4). The relative prey mass of herring in the cod stomachs showed low spatiotemporal variation compared to the other prey apart from very local hotspots from year to year (Figure 3a,d). Both the relative prey mass of and predation on Saduria were highest in the central parts of the southern Baltic Sea (Figures 3b,e and 4b,e), which corresponds to the core area of the Saduria distribution in this region (Appendix S1: Figure S8b). Both the relative prey mass of sprat (Figure 3c,f) and the predation on sprat (Figure 4c,f) occurred throughout the Baltic Sea, although the predation was more limited to the southwestern part in recent years due to the shift in distribution of cod (Figure 5).
Relative prey mass for a cod of 33 cm. Colors indicate the median value across 500 simulated spatial predictions for herring (a, d), Saduria (b, e), and sprat (c, f), for years 1994 (top row, a–c) and 2022 (bottom row, d–f), as examples. The color scale is square‐root‐transformed to better visualize the spatial patterns. Note that scales are shared within species, across years. Only grid cells with depth <130 m are included in the plot.
Cod density × relative prey mass (i.e., predation) plotted in space. Colors indicate the mean across 500 simulated spatial predictions of both relative prey mass and cod density for herring (a, d), Saduria (b, e), and sprat (c, f), for years 1994 (top row, a–c) and 2022 (bottom row, d–f), as examples. The sum of this metric (expanded by grid cell area) across space is the relative population‐level predation intensity depicted in Figure 6. Note that the color scale is 3rd‐root‐transformed and that scales are shared within species, across years. Only grid cells with depth <130 m are included in the plot.
Cod biomass density in space. Colors indicate the mean across 500 simulated spatial predictions for years 1994 (a) and 2022 (b), as examples. Note that the color scale is square‐root‐transformed and truncated at the 99.9th percentile to better visualize the spatial patterns (maximum cod biomass density is 3608 kg/km2). Only grid cells with depth <130 m are included in the plot.
Over time, the area‐expanded per‐capita predation on herring increased steadily (with a small decline around the mid 2000s), while the population‐level predation reached a peak around 2010 after which it declined to levels similar to the early 2000s (Figure 6a,d). Both per‐capita and population‐level predation on Saduria showed a peak at around 2007, after which it declined steadily to very low levels (Figure 6b,e). Both per‐capita and population‐level predation on sprat declined in the early 1990s (Figure 6c,f). From then, the per‐capita predation on sprat varied somewhat cyclically and showed a weak tendency for an overall increase throughout the time period (Figure 6c). The population‐level predation on sprat instead peaked around 2010, and has declined since (Figure 6f). The uncertainty around the per‐capita and total predation estimates is substantial when accounting for uncertainty in both cod biomass density and relative prey masses, especially for Saduria early in the time series and sprat (Figure 6b,c), and this is largely due to higher uncertainty in the stomach content predictions (Appendix S1: Figure S9). Time series of per‐capita and population‐level predation were nearly identical for Saduria and sprat when empty stomachs were omitted. For herring, trends over time tended to be flattened when empty stomachs were omitted (Appendix S1: Figure S10).
Relative per‐capita (top row, a–c) and population‐level predation (bottom row, d–f) by cod on herring (a, d), Saduria (b, e), and sprat (c, f) over time. Points depict the median predation, and vertical lines depict the range between the 10th and 90th percentile of predation, calculated from 500 simulated spatial predictions of both relative prey mass and cod density. Blue lines depict fits from a generalized additive model with year modeled as a penalized spline, and ribbons correspond to the 50% and 90% credible interval of the prediction.
The overlap between cod and its prey was highest in the central parts of the southern Baltic Sea, and along the southeast coast of Sweden (Figure 7). Over time, the overlap with herring started to decline around 2005 (Figure 8a). The spatial overlap with Saduria increased slightly between 1993 and 2009, but since 2010, it has been lower than average, resulting in a weakly negative trend over time (Figure 8b). The spatial overlap with sprat also declined over time, but in a more cyclic fashion. The current spatial overlap with sprat is the lowest since 2007 (Figure 8c).
Spatial overlap between cod and its prey herring (a), Saduria (b), and sprat (c). Colors indicate the mean across 500 simulated spatial predictions of cod density in years 1994 (top row, a–c) and 2022 (bottom row, d–f), as examples. The color scale is 3rd‐root‐transformed to better visualize the spatial patterns. Only grid cells with depth <130 m are included in the plot.
Spatial overlap between cod and its prey herring (a), Saduria (b), and sprat (c). Points depict the median overlap, and vertical lines depict the range between the 10th and 90th percentile of overlap, calculated from 500 simulated spatial predictions of cod density. Blue lines depict fits from a generalized additive model with year modeled as a penalized spline, and ribbons correspond to the 50% and 90% credible interval of the prediction.
The correlation between annual estimates of per‐capita and population‐level predation with spatial overlap was only clearly positive for Saduria (Figure 9b,e). For herring and sprat, the correlation coefficients ranged between −0.39 and −0.12, but the CIs overlapped zero in all cases but the per‐capita predation on herring (Figure 9a,c,d,f). The results were nearly identical when fitting diet models to data with empty stomachs omitted for Saduria and sprat, whereas for herring, the main difference was that the CIs of the correlation between per‐capita herring predation and cod‐herring overlap crossed 0 (Appendix S1: Figure S11).
Correlation between relative per‐capita predation and spatial overlap (top row, a–c), and relative population‐level predation and spatial overlap (bottom row, d–f). Points depict the median, and vertical and horizontal lines depict the range between the 10th and 90th percentile of predation and spatial overlap, respectively, calculated from 500 simulated spatial predictions of predation and cod density. The Pearson correlation coefficient and its 95% CI are printed in the top right corner of each panel for herring (a, d), Saduria (b, e), and sprat (c, f).
DISCUSSION
In this study, we quantified the effects of local prey availability on the stomach contents of predators and the relationship between predation and spatial overlap of predator and prey. Our analysis used the Baltic Sea as a case study. However, the spatiotemporal modeling approach used, which scales up local‐scale diet data to population‐level predation metrics while accounting for range shifts, can be applied generally, and is well suited to improve our understanding of predator–prey interactions in other systems. Our results suggest that the effects of local prey availability are modest and uncertain and that only in the case of the benthic isopod Saduria is availability related to what is found in the stomachs of cod. Furthermore, only for Saduria do we find a correlation between predator–prey spatial overlap and per‐capita and population‐level predation. Our analysis therefore echoes previous studies on the challenges of linking prey dynamics to predator performance in marine generalist predators (Fall et al., 2021; Goodman et al., 2022; Hilborn et al., 2017) and is at contrast with the strong bottom‐up effects present in trophic ecosystem models (Chagaris et al., 2020; Smith et al., 2011). Our analyses constitute important steps toward understanding the spatial scale of trophic interactions (Amarasekare, 2008; Carroll et al., 2024), which are needed to support the implementation of ecosystem‐based fisheries management.
Effects of prey availability on feeding of generalist predators
It has been suggested that the large abundance fluctuations of many small pelagic fish species would help identify relationships between forage fish and predators (Hilborn et al., 2017) because it results in large contrasts for analysis. This is also the case in our system. The Baltic Sea sprat stock increased fivefold between 1991 and 1997 (ICES, 2022b), albeit in the entire Baltic Sea, and not necessarily in the entire distributional range of cod. However, we did not detect an effect of sprat availability on the relative mass of sprat in cod stomachs, and the per‐capita and population‐level predation on sprat by cod correlated negatively (but not significantly) with cod–sprat spatial overlap. One potential reason for this could be that cod are not limited by sprat biomass but by their digestive capacity. This could limit predation even if the population‐level prey biomass is low, as long as there is relatively high prey density locally (Fall & Fiksen, 2020; Hilborn et al., 2017). Alternatively, it may be that the prey availability is not modeled at the appropriate spatial or temporal scale. Prey biomass is aggregated over ICES rectangles and is derived from a hydroacoustic survey conducted in the fourth quarter (while approximately 60% of stomach data are sampled in the fourth quarter). Hence, the spatiotemporal mismatch between data might be too large to accurately reflect predator–prey interactions, given the patchy distribution of schooling prey and the fact that stomach content data reflect consumption over a relatively short time period. Moreover, we measure overlap and pelagic prey availability in two dimensions, which may not accurately reflect encounter rates in a three‐dimensional environment. In line with this, Fall et al. (2021), also found weak effects of prey (capelin, Mallotus villosus) abundance on the consumption of capelin by cod. Instead they found that the proximity of capelin to the seafloor was a better predictor of capelin in the cod diet, which illustrates the potential importance of the third dimension (i.e., depth). Future studies could explore whether this is also the case in the Baltic Sea, potentially using high‐resolution data on schooling fish in combination with stomach content data to try and determine appropriate scales for analysis.
The spatial and temporal scale is not only relevant for the actual interaction (encounter and predation), but also when it comes to “scaling‐up” functional responses or predation metrics from experimental or local scales to scales more relevant for management, for example, to the stock or population level (Hunsicker et al., 2011). This is inherently difficult, as the relationship between consumption and prey density can change over spatial scales (Bergström & Englund, 2004). Another approach, which overcomes this issue, is to use stock‐level estimates when estimating functional responses (Essington & Hansson, 2004). However, a limitation is that spatial heterogeneity in feeding dynamics or predator distribution and range shifts would not be accounted for. In the alternative model‐based approach used here, we arrive at population‐level predation metrics by summing spatially explicit predictions over the heterogeneous domain. This has the benefit that local‐scale heterogeneity in prey availability and stomach contents is explicitly considered, providing more accurate estimates.
Implications for Baltic Sea cod
Our results on the spatiotemporal dynamics of predation and overlap provide novel insights that both corroborate and contrast previous hypotheses on Baltic cod. Neuenfeldt et al. (2020) identified that cod feeding rates on Saduria and sprat were substantially higher in the period 1963–1989 than in 1994–2014. Within the second time period, the physiological condition of cod declined rapidly (Eero et al., 2023; Lindmark et al., 2023; Mion et al., 2021). This makes it an important time period to analyze for understanding the relationship between feeding dynamics and predator–prey overlap. We observe that feeding on Saduria continued to decline after 2014 to near zero in the most recent years. Although the spatial overlap with Saduria also declined slightly, it seems that the loss of predation on Saduria is mainly driven by changes in per‐capita predation, because trends are similar for per‐capita (where the total cod biomass is factored out) and population‐level predation rates. This could in turn be due to declining local abundances of Saduria, because even though Saduria does not seem to have declined substantially over the last 30 years in shallow areas (Svedäng et al., 2022), its area occupied may have declined, since its depth distribution is linked to oxygen dynamics on the sea floor (Karlson et al., 2002). This, however, is not currently possible to investigate further due to the low spatiotemporal resolution of Saduria biomass data. Increased competition for Saduria with flounder (Platichthys spp.) may potentially explain the declines in the availability of Saduria to cod, as hypothesized in Haase et al. (2020). Support for this hypothesis is found in recent years, as high flounder density is associated with lower levels of Saduria in predator stomachs (Lindmark et al., 2023). However, whether competition explains the long‐term decline in Saduria in cod stomachs is not as clear, since in the 1980s and early 1990s (Orio et al., 2017), flounder was more numerous than now and cod still fed on Saduria in high numbers.
Previous studies have hypothesized that the decline in feeding rates on sprat is linked to availability and spatial overlap (Casini et al., 2016; Eero et al., 2012; Neuenfeldt et al., 2020), although this has never been tested explicitly. While there is a spatial mismatch in the sense that sprat biomass is highest in the northeast (Appendix S1: Figure S8c) and cod is mainly found in the southern parts, that does not necessarily mean cod is limited by low sprat abundances in the south. In line with Eero et al. (2012), our results suggest that predation pressure on sprat is highest in the south. We find fluctuations in the spatial overlap with sprat, with a positive trend for per‐capita predation. Also considering the lack of effect of sprat biomass on the relative mass of sprat in cod stomachs (Table 1), our study does not support the hypothesis that spatial dynamics and spatiotemporal mismatch in overlap explain trends in cod feeding on sprat. The lack of correlation between overlap and predation is also found in herring, and per‐capita predation on herring by cod even increased since 2010 when the spatial overlap started to decline.
Limitations and areas of future research
We made several simplifying assumptions that warrant future research. When calculating the predation indices, we used the total cod density in the survey catches, that is, not resolved by length class. Since cod undergo strong ontogenetic diet shifts (Kulatska et al., 2019; Lindmark et al., 2025), it may be more appropriate to use the biomass density of cod within specific size groups that mainly feed on a specific prey. However, that may also not have a strong impact, since the correlation between local cod densities of different size groups is quite high (Jacobsen et al., 2023), the population size structure is highly truncated, and the main sizes currently in the population feed on a mix of pelagic and benthic prey (Lindmark et al., 2025). Moreover, the predation indices calculated here are only for cod of 33 cm (for simplicity, the mean size in the diet data), although they could be predicted for any size since length is a covariate in the models. This simplifying assumption effectively assumes the entire cod stock is of a specific size, which can be misleading for population‐level predation metrics since the size distribution of the cod population has been truncated over time (Eero et al., 2023). However, it does facilitate comparison over time for the per‐capita predation metrics. Moreover, the indices are relative since the biomass density of cod is relative (due to unknown survey catchability), and since length is only a fixed effect in the model, predicting for different lengths would just change the relative value of the index identically over time.
The exact values of the predation indices are not directly comparable to mortality rates because the units are different. However, it would be straightforward to expand our work by converting our predation metrics to predation rates using gastric evacuation models (as in e.g., Gartland & Latour, 2024; Tengvall et al., 2024). A qualitative comparison between our predation estimates and predation estimates from the Baltic Stochastic Multi‐Species model (SMS) (Lewy & Vinther, 2004) suggests that population‐level predation intensity and natural mortality of sprat and herring both increased from 1993 until around 2010, after which cod predation declined to levels comparable to around the year 2000 (ICES, 2019). This implies that the stomach content data, more than model type (spatiotemporal index or multispecies assessment model) drives the result, although more work is needed to identify which processes cause discrepancies between the time series of predation and mortality.
Stomach content sampling is often affected by gaps and inconsistencies in space and time (Figure 1). This represents a challenge for using stomach data in population dynamic models without some form of standardization, since there are clear spatiotemporal patterns in stomach content data. This is, however, not routinely done (ICES, 2019; Neuenfeldt et al., 2020). The approach used in this study represents a model‐based approach to estimate trophodynamic indices over the spatial domain using spatial and spatiotemporal random effects (Cao et al., 2017; Karp et al., 2025; Thorson et al., 2015). This approach has benefits over design‐based indices by being able to include covariates and latent variables, which can improve estimates when the sampling is spatially or temporally unbalanced as the Baltic diet data are (Figure 1). This means estimates in areas with low sampling intensity for a given year depend on the ability to estimate constant and time‐varying spatial random effects, which illustrates the importance of critically evaluating models (Yalcin et al., 2023).
We believe there are several possibilities for future work. For instance, this could include further exploration of density dependence, or at which spatiotemporal scale covariates should be included (Fall et al., 2021; Lindmark et al., 2025). Another area of research could be the covariation between different prey groups, since they are likely not independent. For instance, if cod recently fed on herring, they may not feed on benthic prey soon after, or cod may switch to sprat if the abundance of herring declines. Questions such as these could potentially be addressed using dynamic structural equation models including simultaneous or lagged effects (Thorson et al., 2024).
CONCLUSIONS
Predator–prey interactions play an important role in defining ecosystem functioning and the trophic structure of marine food webs. The relationship between predator feeding dynamics and prey is a crucial aspect of this. Understanding these interactions is critical for supporting the implementation of broader food web considerations in an ecosystem‐based approach to fisheries management. Essentially, this means being able to predict the ecological effects on predators that may stem from changes in prey abundance and distribution, fisheries and other anthropogenic pressures. For increases in prey to have a causal impact on predator productivity, there has to be a link between prey availability and predation. While our analysis does not control for all potential confounding factors, we do not find evidence of such positive associations between pelagic prey and predator feeding dynamics despite large fluctuations in abundance, local‐scale availability, and population‐level spatiotemporal overlap. While we do acknowledge that spatiotemporal dynamics of these interactions are complex and scale‐dependent and may be difficult to quantify (e.g., Fall et al., 2021), our results could mean that the effects of specific prey are weaker than previously thought. However, it is difficult to make general statements given the mixed results in the literature (e.g., Free et al., 2021; Goodman et al., 2022), and since it may have implications for assessment and management, it should be evaluated case by case. In the example of Baltic cod, it could mean that management interventions aimed at increasing the availability of pelagic prey would have limited impacts on the productivity of cod.
AUTHOR CONTRIBUTIONS
Conceptualization: Max Lindmark, Christopher A. Griffiths, Valerio Bartolino, and Federico Maioli. Data curation: Viktor Thunell, Max Lindmark. Methodology: Max Lindmark, Christopher A. Griffiths, Valerio Bartolino, Federico Maioli, and Sean C. Anderson. Formal analysis: Max Lindmark, Viktor Thunell, and Sean C. Anderson. Visualization: Max Lindmark and Sean C. Anderson. Software: Sean C. Anderson. Project administration: Nis S. Jacobsen, Max Lindmark, Christopher A. Griffiths, and Valerio Bartolino. Resources: Katarzyna Nadolna‐Ałtyn, Joanna Pawlak, Marzenna Pachur, Marcin Rakowski, Karolina Wikström, Mayya Gogina, and Didzis Ustups. Writing—original draft: Max Lindmark. Writing—review and editing: All authors.
CONFLICT OF INTEREST STATEMENT
The authors declare no conflicts of interest.
Supporting information
Appendix S1.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Aitchison, J. 1955. “On the Distribution of a Positive Random Variable Having a Discrete Probability Mass at the Origin.” Journal of the American Statistical Association 50: 901. 10.2307/2281175 · doi ↗
- 2Akaike, H. 1973. “Information Theory and an Extension of the Maximum Likelihood Principle.” In Second International Symposium on Information Theory, edited by B. N. Petrov and F. Caski , 267–281. Budapest: Akademiai Kiado.
- 3Amarasekare, P. 2008. “Spatial Dynamics of Foodwebs.” Annual Review of Ecology, Evolution, and Systematics 39: 479–500. 10.1146/annurev.ecolsys.39.110707.173434. · doi ↗
- 4Anderson, S. C. , E. J. Ward , P. A. English , L. A. K. Barnett , and J. T. Thorson . 2024. “sdm TMB: An R Package for Fast, Flexible, and User‐Friendly Generalized Linear Mixed Effects Models with Spatial and Spatiotemporal Random Fields.” bio Rxiv 2022.03.24.485545. 10.1101/2022.03.24.485545. · doi ↗
- 5Bakun, A. 2006. “Wasp‐Waist Populations and Marine Ecosystem Dynamics: Navigating the “Predator Pit” Topographies.” Progress in Oceanography 68: 271–288. 10.1016/j.pocean.2006.02.004. · doi ↗
- 6Bartolino, V. , H. Tian , U. Bergström , P. Jounela , E. Aro , C. Dieterich , H. E. M. Meier , M. Cardinale , B. Bland , and M. Casini . 2017. “Spatio‐Temporal Dynamics of a Fish Predator: Density‐Dependent and Hydrographic Effects on Baltic Sea Cod Population.” P Lo S One 12: e 0172004. 10.1371/journal.pone.0172004.28207804 PMC 5313222 · doi ↗ · pubmed ↗
- 7Bergström, U. , and G. Englund . 2004. “Spatial Scale, Heterogeneity and Functional Responses.” Journal of Animal Ecology 73: 487–493. 10.1111/j.0021-8790.2004.00823.x. · doi ↗
- 8Bürkner, P. C. 2017. “brms: An R Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software 80(1):1–28. 10.18637/jss.v 080.i 01 · doi ↗