2022 Bluefish (Pomatomus saltatrix) Research Track Working Paper 4: Forage Index

Vector Autoregressive Spatio-Temporal (VAST) modeling of piscivore stomach contents, 1985-2021

Abstract

Does prey drive availability of bluefish? Assessing small pelagic fish trends in space and time using piscivore diet data

Changing distribution and abundance of small pelagics may drive changes in predator distributions, affecting predator availability to fisheries and surveys. However, small pelagic fish are difficult to survey directly, so we developed a novel method of assessing small pelagic fish aggregate abundance via predator diet data. We used piscivore diet data collected from multiple bottom trawl surveys within a Vector Autoregressive Spatio-Temporal (VAST) model to assess trends of small pelagics on the Northeast US shelf. The goal was to develop a spatial “forage index” to inform survey and/or fishery availability in the bluefish (Pomatomus saltatrix) stock assessment. Using spring and fall surveys from 1973-2020, 20 small pelagic groups were identified as major bluefish prey using the diet data. Then, predators were grouped by diet similarity to identify 21 piscivore species with the most similar diet to bluefish in the region. Diets from all 22 piscivores were combined for the 20 prey groups at each surveyed location, and the total weight of small pelagic prey per predator stomach at each location was input into a Poisson-link delta model to estimate expected prey mass per predator stomach. Best fit models included spatial and spatio-temporal random effects, with predator mean length, number of predator species, and sea surface temperature as catchability covariates. Spring and fall prey indices were split into inshore and offshore areas to reflect changing prey availability over time in areas available to the recreational fishery and the bottom trawl survey, and also to contribute to regional ecosystem reporting.

1 Introduction

The objective of this work was to create a “prey index” to evaluate changes in bluefish prey over time and in space using Vector Autoregressive Spatio-Temporal modeling (VAST (Thorson and Barnett 2017; Thorson 2019)). This approach was patterned on Ng et al. (2021), which used predator stomach data to create a biomass index for a single prey, Atlantic herring. Expected biomass of herring per stomach was estimated in VAST with 2 linear predictors, the number of herring per stomach and the average weight of herring in a stomach.

We adapted the approach of Ng et al. (2021) to get an index for “bluefish prey” in aggregate rather than a single prey species. Further, we include inshore and offshore regions by combining results across regional bottom trawl surveys surveys as was done for summer flounder biomass in Perretti and Thorson (2019). Finally, since bluefish themselves are somewhat sparsely sampled by the surveys, we aggregate all predators that have a similar diet composition to bluefish to better represent bluefish prey biomass.

We characterize mean weight of bluefish prey from all piscivores caught at each survey location and model that over time/space. Covariates potentially affecting perceived patterns in the bluefish prey index include number of predators, size composition of predators, and sea surface temperature (SST) at each survey location.

Therefore, the steps involved to estimate the forage index included defining the input dataset, and running multiple configurations of the VAST model. Steps involved in defining the dataset included defining “bluefish prey”, defining a set of piscivore predators with similar diets to bluefish, integrating diet data from two regional surveys, and integrating supplementary SST data to fill gaps in in-situ temperature data measurements. Steps involved in running the VAST model included decisions on spatial footprint, model structure, model selection to determine if spatial and spatio-temporal random effects were supported by the data, and further model selection to determine which catchability covariates were best supported by the data. Finally, subsets of the spatial domain were defined to match bluefish assessment inputs (survey and recreational fishery CPUE) for potential use as covariates in bluefish stock assessment models, and a bias-corrected (Thorson and Kristensen 2016) forage index for each spatial subset was generated.

This approach is generalizable to any spatial subset of the full VAST spatial domain, such that forage indices for each ecoregion on the Northeast US continental shelf, or for other piscivore predators can also be generated.

2 Methods

2.1 Input dataset

Fish food habits data are collected aboard several regional fishery independent surveys in the Northeast US. The longest time series of diet has been collected by the Northeast Fisheries Science Center (NEFSC) bottom trawl survey (Smith and Link 2010) from south of Cape Hatteras, NC to the Scotian Shelf at the US/Canada border since the early 1970s, which represents the bulk of the data for this analysis. Using similar survey protocols to the NEFSC survey, the NorthEast Area Monitoring and Assessment Program (NEAMAP) survey has collected fish diet data in inshore waters from Cape Hatteras, NC to Cape Cod, MA since 2008 (Northeast Area Monitoring & Assessment Program (NEAMAP) et al. 2021). All analyses were completed in R (Team 2021).

2.1.1 Defining bluefish prey

Bluefish eat small pelagics that are not well sampled by bottom trawl surveys. Bluefish themselves are not well sampled by bottom trawl surveys. Nevertheless, the diet samples collected for bluefish indicate that anchovies, herrings, squids, butterfish, scup, and small hakes are important prey. See the Bluefish Ecosystem Socioeconomic Profile working paper, as well as this web summary page link for NEFSC survey and this presentation link for NEAMAP survey bluefish diet composition summaries.

Using all sampled bluefish stomachs included in the NEFSC food habits database 1973-2021, we created a list of all pelagic nekton bluefish prey that had at least 10 observations (Table 2.1). Broad categories such as empty stomach, fish unidentified, Osteichthyes, unidentified animal remains, and blown stomach were not included in the prey list.

2.1.2 Defining piscivore predators

The choice of predators is largely intended to balance increasing sample size for modeling bluefish prey with using predators likely to be foraging similarly to bluefish. One extreme assumption would be to include only bluefish as predators, but there are relatively few bluefish diet samples due to incomplete availability to bottom trawl surveys (see supplemental information). This would miss prey available to bluefish because we have not sampled bluefish adequately. The opposite extreme assumption would be to include all stomachs that contain any of the top bluefish prey, regardless of which species ate the prey. This would include predators that do not forage similarly to bluefish and might therefore “count” prey that are not actually available to bluefish due to habitat differences. The intermediate approach which selects a group of piscivores that forage similarly to bluefish both increases sample size and screens out the most dissimilar predators to bluefish. A further refinement to the input data is only using the prey items identified above as “bluefish prey” across all predators identified as piscivores.

For bluefish forage index modeling, we are selecting a set of predators that have high diet similarity to bluefish. Garrison and Link (2000) evaluated similarity of predator diets on the Northeast US shelf to develop foraging guilds. The Schoener similarity index (Schoener 1970) was applied

to assess the dietary overlap, \(D_{ij}\), between predator pairs:

\[ D_{i,j} = 1 – 0.5 (\sum |p_{i,k} – p_{j,k}|) \]

where \(p_{i,k}\) = mean proportional volume of prey type \(k\) in predator \(i\) and \(p_{i,k}\) = mean proportional volume of prey type \(k\) in predator \(j\) (Garrison and Link 2000).

Garrison and Link (2000) used NEFSC bottom trawl survey data 1973-1997. We are using diets from 1985-2020 to characterize the forage index. Therefore an additional 20+ years of diet information is available to assess whether predator diet similarity has changed. Diet similarity analysis has been completed (B. Smith, pers. comm.), with a table of prey similarity available via this link on the NEFSC food habits shiny app that was used to define feeding guilds based on 50 predators. We evaluated results from several clustering algorithms (see this link) to determine how robustly different sets of predators grouped with bluefish. The working group selected the piscivore list based on the “complete” clustering algorithm, including all species that clustered with all 3 sizes of bluefish (Table 2.2).

This piscivore dataset better captured predators that feed similarly to bluefish (e.g. striped bass), and has a higher proportion of stations with bluefish prey than a dataset based on the Garrison and Link (2000) piscivore definition. We also evaluated a piscivore definition using only the predators that always cluster with bluefish no matter what clustering algorithm is applied. However, a dataset based on that limited piscivore list excluded predators highlighted by bluefish experts (e.g., striped bass) and resulted in the inclusion of fewer overall stations than the either of the above piscivore definitions, with a lower proportion of included stations with bluefish prey.

The NEAMAP survey operates closer to shore than the current NEFSC survey. While both surveys capture many of the same predators, some are not available close to shore and others are more available close to shore. The NEAMAP dataset includes the following piscivore predators and size ranges, adding two species not captured by the NEFSC survey offshore but included based on working group expert judgement of prey similarity to bluefish (Spanish mackerel and spotted sea trout) and leaving out those not captured inshore:

• Summer Flounder 21-70 cm
• Silver Hake 21-76 cm
• Weakfish 26-50 cm
• Atlantic Cod 81-150 cm
• Bluefish 3 – 118 cm
• Striped Bass 31 – 120 cm
• Spanish Mackerel 10 – 33.5 cm
• Spotted Sea Trout 15.5 – 34 cm
• Spiny Dogfish 36 – 117 cm
• Goosefish 5 – 124 cm

2.1.3 Integrating regional surveys

For each survey dataset, the full diet database was filtered to include only predators on the list of piscivores with the most diet similarity to bluefish (Table 2.2). Then, the list of bluefish prey (Table 2.1) was used to categorize prey items for each predator as “bluefish prey” or “other prey”. Each station was given a unique station identifier (cruise and station number), and the total weight (g) of bluefish prey at each station was summed. Total bluefish prey weight was divided by the total number of stomachs across all piscivore predators at the station to get mean bluefish prey weight (g) at each station. In addition, the number of piscivore species and the mean size (total length, cm) across all piscivores was calculated at each station. Seasons were identified as “Spring” (collection months January - June) and “Fall” (collection months July-December) to align with assumptions used in the bluefish stock assessment. Vessel identifiers were assigned based on years and survey, with two vessels used for the NEFSC survey (R/V Albatross prior to 2009 and R/V Bigelow 2009 to present) and a single vessel used for the NEAMAP survey 2008-present. These were used to evaluate whether vessel effects were present. Variable names were reconciled between NEFSC and NEAMAP, and the datasets were appended into a single dataset with one row per station including station ID, year, season, date, latitude, longitude, vessel, mean bluefish prey weight, mean piscivore length, number of piscivore species, and sea surface temperature (degrees C).

2.1.4 Filling gaps in sea surface temperature (SST) data

Initial dataset checks revealed that approximately 10% of survey stations were missing in-situ sea water temperature measurements. Gaps in temperature information were more prevalent early in the time series (1980s and early 1990s), although stations without temperature data were found in nearly all years (see “SST mismatch by year” section at this link). Rather than truncate the dataset to only those stations with in-situ temperature measurements, we investigated other sources of SST data to fill gaps.

Two SST data sources were investigated, both based on satellite data: the National Oceanic and Atmospheric Administration Optimum Interpolation Sea Surface Temperature (NOAA OI SST) V2 High Resolution Dataset (Reynolds et al. 2007) data provided by the NOAA PSL, Boulder, Colorado, USA, from their website at https://psl.noaa.gov, and the higher resolution source data for the OI SST, the AVHRR Pathfinder SST data linked here (Saha et al. 2018). Both sources provide global daily SST at different spatial resolutions (OI SST uses a 25 km grid, and AVHRR uses a 4 km grid) from 1981-present.

The OI SST data are provided as global files for each year. Files for years 1985-2021 were downloaded from https://downloads.psl.noaa.gov/Datasets/noaa.oisst.v2.highres/sst.day.mean.[year].v2.nc as rasters using code developed by Kim Bastille for Northeast US ecosystem reporting, cropped to the Northeast US spatial extent, and converted to R dataframe objects where the temperature of a grid cell is associated with the coordinates at the center of the grid cell. Then, OI SST temperature was matched to the survey data using year, month, day and spatial nearest neighbor matches to survey station locations.

AVHRR data is organized into year/data folders with 2 nc files for each date, one day and one night. A list of survey year-month-day combinations was generated from the piscivore diet dataset to download only AVHRR daily files within \(\pm\) 2 days of survey dates from https://www.ncei.noaa.gov/data/oceans/pathfinder/Version5.3/L3C/[year]/data/. The “sea_surface_temperature” and “quality_level” fields were extracted as rasters, cropped to the Northeast US spatial extent, appended into an annual dataset and saved as R dataframe objects similar to the OI SST methods. Examination of a subset of AVHRR SST data with quality level 3 and above showed extended periods with little SST data within the survey domain, likely due to cloud cover (see “SST AVHRR 2021 survey dates” section at this link). Therefore, to match AVHRR SST to survey stations, methods for combining and filling temporal and spatial gaps would be required, which was beyond the scope of the current analysis.

For survey stations with in-situ temperature measurements, the in-situ measurement was retained. For survey stations with missing temperature data (~10% of all stations), OI SST was substituted for input into VAST models.

2.2 VAST modeling

We used VAST (Thorson and Barnett 2017; Thorson 2019) to evaluate changes in bluefish prey biomass and distribution over time. VAST is structured to estimate fixed and random effects across two linear predictors, which are then multiplied to estimate an index of the quantity of interest. Using notation from Thorson (2019), a full model for the first linear predictor \(\rho_1\) for each observation (\(i\)) can include fixed intercepts (\(\beta\)) for each category (\(c\)) and time (\(t\)), spatial random effects (\(\omega\)) for each location (\(s\)) and time, spatio-temporal random effects (\(\varepsilon\)) for each location, category, and time, fixed vessel effects (\(\eta\)) by vessel (\(v\)) and category, and fixed catchability impacts (\(\lambda\)) of covariates (\(Q\)) for each observation and variable (\(k\)):

\[ \rho_1(i) = \beta_1(c_i, t_i) + \omega_1^*(s_i, c_i) + \varepsilon_1^*(s_i, c_i, t_i) + \eta_1(v_i, c_i) + \sum_{k=1}^{n_k} \lambda_1(k) Q(i,k)\]

The full model for the second linear predictor \(\rho_2\) has the same structure, estimating \(\beta_2\), \(\omega_2\), \(\varepsilon_2\), \(\eta_2\), and \(\lambda_2\) using the observations, categories, locations, times, and covariates. VAST models can also include habitat (density) covariates, which we did not implement here, and have left out of the equation for simplicity.

2.2.1 Structural decisions

Thorson (2019) lists 15 major decisions for constructing a VAST model. These include decisions on spatial domain, categories modeled (species, ages, etc), data type (presence absence, number, weight), including spatial and or spatio-temporal variation, spatial resolution, univariate vs multivariate response and factors, specifying temporal correlation, including density and or catchability covariates, treatment of area swept calculation, including vessel effects, selecting a link function, specifying derived quantities, bias correcting derived quantities, and model selection. Here we outline the decisions made in developing the forage index.

2.2.1.1 Spatial domain

Models were run using the full Northwest Atlantic grid built into VAST (Fig. 2.1). Specific strata sets were used from this full model to develop indices matching the spatial extent of different assessment inputs (see below).

Figure 2.1: Northwest Atlantic Grid (blue outline) from https://github.com/James-Thorson-NOAA/FishStatsUtils

2.2.1.2 Categories, data type, and link function

We model all bluefish prey in aggregate as a single category. The mean weight of bluefish prey in a stomach at each location is treated as biomass data in the model. Therefore, VAST applies a delta model where the first linear predictor models encounter rate and the second linear predictor models amount of prey (equivalent to positive catch rates on a survey). Following what Ng et al. (2021) did for herring, as well as recommended practice for biomass data (Thorson 2019), we apply a Poisson-link delta model to estimate expected prey mass per predator stomach.

2.2.1.3 Spatial variation, resolution, response type, and temporal correlation

We include spatial and spatio-temporal variation in both linear predictors, but test whether the data support these effects using model selection (see below). Similar to Ng et al. (2021) we used the default spatial smoother in VAST, the stochastic partial differential equation (SPDE) approximation to the Mat'ern correlation function (method = “mesh”; Thorson (2019)). Although directional correlation (anisotopy) can be common in fishery collections with depth gradients along a continental shelf (Thorson 2019), we tested whether the inclusion of anisotopy as a fixed effect was supported using model selection (see below). We used a higher spatial resolution than Ng et al. (2021) did for herring-predator pairs, here defining 500 “knots” or standardized locations optimally allocated among all observed survey stations in the full dataset as estimated by k-means clustering of the data, to define the spatial dimensions of each seasonal model and the annual model.

Modeling all bluefish prey in aggregate leads to a univariate model, producing a single forage index which is most easily integrated into the bluefish assessment model. We did not include temporal correlation in fixed intercepts to maintain independence of forage abundance in each modeled year (Thorson 2019). We did not include temporal correlation in spatio-temporal random effects because most survey areas were sampled each year, so projecting forage “hotspots” between years using temporal correlation was not desirable for this application.

2.2.1.4 Including covariates, vessel effects, area swept, and other decisions

We explored multiple combinations of catchability covariates and vessel effects. Surveys were conducted aboard multiple vessels over time and between NEFSC and NEAMAP, so we investigated vessel effects for the NEAMAP vessel and NEFSC vessels Albatross and Bigelow (commonly included in regional stock assessments when survey indices are not modeled separately). Vessel effects were modeled as overdispersion parameters. Catchability covariates explored included mean predator length at each station, number of predator species at each station, and sea surface temperature (SST) at each station.

The predator length covariate may more directly model vessel differences in predator survey catch that affect stomach contents than modeling a vessel catchability covariate directly. Ng et al. (2021) found that predator length covariates were strongly supported as catchability covariates (larger predators being more likely to have more prey in stomachs). The rationale for including number of predator species is that more species “sampling” the prey field at a particular station may result in a higher encounter rate (more stations with positive bluefish prey in stomachs). Water temperature was also included as a potential catchability covariate, because temperature affects predator feeding rate.

VAST can include habitat or density covariates that are expected to drive modeled species distribution and abundance (as opposed to catchability covariates, which affect our survey observations). For example, a certain habitat or depth may limit the range or productivity of a species. Because we are interested in an aggregate index of forage fish that includes a diversity of species that use many different habitats, including density covariates appropriate across all species (that affect density in the same way) may not be feasible, and was not explored for this project.

Although the forage fish index is based on trawl-surveyed fish predators and the area swept of the net capturing predators is available, determining the actual area swept of the predators “sampling” the prey field is less clear. Therefore, we set area swept to 1 as recommended for “sampling gears” with unknown effective sampling areas, which means our forage abundance index does not have an interpretable scale, but should be proportional to actual forage biomass (Thorson 2019).

The derived quantity of interest here is a biomass index for each of spring, fall, and annual datasets for bluefish prey species. We have also included supplementary plots of the center of gravity for each seasonal model and the annual model.

Bias correction of the forage fish biomass index for each model (and spatial subdivisions, see below) is based on Thorson and Kristensen (2016), as implemented in the VAST development branch code (https://github.com/James-Thorson-NOAA/VAST/tree/dev).

2.2.2 Model selection

We compared the AIC for the 500 knot models to see if including the spatial and spatio-temporal random effects in the first and second linear predictors improved the model fit. Model structures tested include with and without anisotopy (2 fixed parameters), and with and without spatial and spatio-temporal random effects in the second linear predictor or both linear predictors. This follows the model selection process outlined in Ng et al. (2021) using restricted maximum likelihood (REML; Zuur et al. (2009)).

Following this, we evaluated catchability covariates using AIC to determine which are best supported by the data. We used the structure selected by the REML model selection, then evaluated vessel effects (overdispersion) and a range of potential catchability covariates using maximum likelihood to calculate AIC instead of REML, because the spatial and spatio-temporal random effects are the same across models.

Our two-step model selection (1: spatial and spatio-temporal random effects, 2: catchability covariates) was completed using the script bluefishdiet/VASTunivariate_bfp_modselection.R.

2.3 Spatial definitions

Our main goal is to determine whether bluefish prey availability has changed in inshore waters where the recreational fishery primarily operates. Our food habits datasets do not extend into inland waters such as bays and sounds, with the exception of Cape Cod Bay. (In the future we might be able to investigate food habits from ChesMMAP or surveys south of Cape Hatteras.) However, there is data from both historical NEFSC surveys and NEAMAP in state coastal waters (0-3 miles from shore), and offshore across the continental shelf.

The model has been partitioned into several definitions of “inshore” and “offshore” for the stock assessment inputs. First we define a partition that is the MAB and GB areas only as the GOM is not relevant to the bluefish assessment. This is called MABGB (Fig. 2.2), while the full model is AllEPU. Within this partition,

1. Survey inshore vs offshore to evaluate availability to the survey index. Strata partitions include:
   • Albatross inshore stations (Fig. 2.3)
   • Bigelow inshore bluefish index stations (Fig. 2.4)
   • offshore bluefish index stations (considered for addition in 2022, Fig. 2.5)
   • offshore non-bluefish stations
2. Recreational fishery inshore vs offshore to evaluate availability to the MRIP CPUE index. Strata partitions include
   • shoreline to 3 miles out (State waters, Fig. 2.6)
   • offshore of 3 miles (Federal waters)

NEFSC survey strata definitions are built into the VAST northwest-atlantic extrapolation grid already. We defined additional new strata to address the recreational inshore-offshore 3 mile boundary. The area within and outside 3 miles of shore was defined using the sf R package as a 3 nautical mile (approximated as 5.556 km) buffer from a high resolution coastline from thernaturalearth R package. This buffer was then intersected with the current FishStatsUtils::northwest_atlantic_grid built into VAST and saved using code here. Then, the new State and Federal waters strata were used to split NEFSC survey strata where applicable, and the new full set of strata were used along with a modified function from FishStatsUtils::Prepare_NWA_Extrapolation_Data_Fn to build a custom extrapolation grid for VAST as described in detail here.

All strata were applied in both seasonal and annual models.

Model output strata of interest for the bluefish assessment are mapped below:

2.3.1 Key Strata

2.3.1.1 Mid-Atlantic and Georges Bank

Figure 2.2: Map of Mid-Atlantic and Georges Bank region.

2.3.1.2 Albatross inshore survey strata

Figure 2.3: Map of inshore survey strata formerly surveyed by the R/V Albatross, now surveyed by NEAMAP.

2.3.1.3 Bluefish Inshore Survey Strata (Bigelow)

Figure 2.4: Map of bluefish inshore survey strata accessible to the R/V Bigelow

2.3.1.4 Bluefish Offshore Survey Strata (proposed 2022)

Figure 2.5: Map of bluefish offshore survey strata proposed for the 2022 assessment

2.3.1.5 State waters

Figure 2.6: Map of state waters: coastline out to 3 nautical miles from shore.

Seasonal models were run using the script bluefishdiet/VASTunivariate_bfp_allsurvs_lencovSST_ALLinoffsplits.R, which contains all stratum definitions. The annual model was run using the script bluefishdiet/VASTunivariate_bfp_allsurvsANNUA:_lencovSST_ALLinoffsplits.R.

The final model runs included all selected covariates, stratum definitions, and bias correction for the biomass index.

3 Results

3.1 Input dataset overview

The list of bluefish prey derived from the most common identifiable prey items in NEFSC diet database (Table 2.1) includes the majority of bluefish diet composition by decade (Fig. 3.1) and season (Fig. 3.2). Colors in the plots show included prey, while gray sections represent “fish unidentified” and other categories not included in this analysis. Even when evaluated annually and seasonally (where bluefish sample sizes may be small for a year-season combination), a majority of observed diet is included in the dataset for analysis (Fig. 3.3).

Figure 3.1: Bluefish diet by decade, NEFSC bottom trawl surveys.

Figure 3.2: Bluefish diet by season, NEFSC bottom trawl surveys.

Figure 3.3: Bluefish diet by season and year, NEFSC bottom trawl surveys.

The full NEFSC diet database 1985-2021 contains 25634 survey stations with diet collections. When including only piscivores feeding similarly to bluefish, the survey stations with diet collections in this time period is 22751. Of these piscivore diet stations, 9027 included our defined bluefish prey, or 39.6773768%. For comparison, 1814 stations have diet samples for bluefish alone, with 905 or 49.8897464% including our defined bluefish prey.

NEAMAP survey stations with diet collections for piscivores (n = 3838) had a higher proportion with our defined bluefish prey (n = 2418, 63.0015633%).

The number of survey stations missing surface temperature data varied considerably by decade. A large percentage of survey stations lacked in-situ temperature measurements between 1985 and 1990, while the percentage of stations missing temperature was generally below 10% (with a few exceptions) from 1991-2021 (Table 3.1). Therefore, OISST temperature estimates were more commonly substituted early in the time series.

3.2 VAST model selection

Comparisons of AIC are presented for both the first (spatial and spatio-temporal random effects, Table ??) and second (catchability covariates, Table ??) rounds of model selection.

Models compared using REML are identified by model name (“modname” in Tables ??) which includes all prey aggregated (“allagg” for all models), season (“all” for annual models of months 1-12, “fall” for models of months 7-12, and “spring” for models of months 1-6), number of knots (500 for all models), and which fixed and random spatial and spatio-temporal effects were included in which linear predictor (1 or 2). The names for model options and associated VAST model settings are:

Model selection 1 (spatial, spatio-temporal effects, no covariates) options (following Ng et al. (2021)) and naming: * All models set Use_REML = TRUE in fit_model specifications.
* Modeled effects, model name suffix, and VAST settings by model:

1. “_alleffectson” = Spatial and spatio-temporal random effects plus anisotropy in both linear predictors; FieldConfig default (all IID)
2. “_noaniso” = Spatial and spatio-temporal random effects in both linear predictors with anisotopy turned off; FieldConfig default (all IID) and use_anisotopy = FALSE
3. “_noomeps2” = Spatial and spatio-temporal random effects plus anisotropy only in linear predictor 1; FieldConfig = 0 for Omega2, Epsilon2
4. “_noomeps2_noaniso” = Spatial and spatio-temporal random effects only in linear predictor 1 with anisotopy turned off; FieldConfig = 0 for Omega2, Epsilon2 and use_anisotopy = FALSE
5. “_noomeps2_noeps1” = Spatial random effects plus anisotropy only in linear predictor 1; FieldConfig = 0 for Omega2, Epsilon2, Epsilon1
6. “_noomeps2_noeps1_noaniso” = Spatial random effects only in linear predictor 1 with anisotopy turned off; FieldConfig = 0 for Omega2, Epsilon2, Epsilon1 and use_anisotopy = FALSE
7. “_noomeps12” = Anisotropy, but no spatial or spatio-temporal random effects in either linear predictor; FieldConfig = 0 for both Omega, Epsilon
8. “_noomeps12_noaniso” = No spatial or spatio-temporal random effects in either linear predictor; FieldConfig = 0 for both Omega, Epsilon and use_anisotopy = FALSE

Using REML, models including spatial and spatio-temporal random effects as well as anisotropy were best supported by the data. This was true for the spring dataset, the fall dataset, and the annual (seasons combined) dataset.

For the second round of model selection with different combinations of vessel effects and or catchability covariates, “modname” in ?? follows a similar pattern as above, with all prey aggregated (“allagg” for all models), season (“all” for annual models of months 1-12, “fall” for models of months 7-12, and “spring” for models of months 1-6), number of knots (500 for all models), and which vessel effects or covariates were included. The names for model options and associated VAST model settings are:

Model selection 2 (covariates) options, FieldConfig default (all IID), with anisotropy:

1. “_base” = No vessel overdispersion or length/number covariates
   
2. “_len” = Predator mean length covariate
3. “_num” = Number of predator species covariate
4. “_lennum” = Predator mean length and number of predator species covariates
5. “_sst” = Combined in situ and OISST covariate
6. “_lensst” = Predator mean length and SST covariates
7. “_numsst” = Number of predator species and SST covariates
8. “_lennumsst” = Predator mean length, number of predator species, and SST covariates
9. “_eta10” = Overdispersion (vessel effect) in first linear predictor (prey encounter)
10. “_eta11” = Overdispersion (vessel effect) in both linear predictors (prey encounter and weight)

Catchability covariates were better supported by the data than vessel effects. Model comparisons above and here led us to the best model fit using mean predator length, number of predator species, and SST at a survey station as catchability covariates.

3.3 Bias-corrected spatial forage indices

3.3.1 Fall Index

Many indices for different spatial strata can be derived from the VAST model. Fall indices (Fig. 3.4) for the full spatial domain (AllEPU) through smaller strata were produced. The same regions are used for all models.

Figure 3.4: All Fall forage indices

We can compare the indices relevant to bluefish on the same scale: inshore (Albatross strata), inshore bluefish, offshore bluefish, and further out, plus the state and federal waters split (Fig. 3.5).

Figure 3.5: Bluefish-related Fall forage indices

Proportions of prey in each bluefish relevant area can be compared (here as a proportion of the full MABGB index; Fig. 3.6).

Figure 3.6: Fall forage indices as proportion of the Mid Atlantic and Georges Bank area.

3.3.2 Fall predicted ln-density

The VAST model predicts density of forage fish across the entire model domain for each year (Fig. 3.7).

Figure 3.7: Yearly maps of VAST model estimated forage fish density for Fall (months 7-12).

3.3.3 Fall Diagnostics

Basic VAST diagnostics include maps of the spatial grid knot placement (“Data_and_knots”), maps of included station locations for each year (“Data_by_year”), residual plots (“quantile residuals”), maps of residuals for each station (“quantile_residuals_on_map”), an anisotropy plot indicating directional correlation (“Aniso”), and a plot of the estimated change in forage fish center of gravity over time (“center_of_gravity”). We present these plots for each of the models in numbered sections below.

3.3.3.1 Data_and_knots

Fall diagnostics:Data_and_knots

3.3.3.2 Data_by_year

Fall diagnostics:Data_by_year

3.3.3.3 quantile_residuals

Fall diagnostics:quantile_residuals

3.3.3.4 quantile_residuals_on_map

Fall diagnostics:quantile_residuals_on_map

3.3.3.5 Aniso

Fall diagnostics:Aniso

3.3.3.6 center_of_gravity

Fall diagnostics:center_of_gravity

3.3.4 Spring Index

Spring indices (Fig. 3.8) for the full spatial domain (AllEPU) through smaller strata were produced.

Figure 3.8: All Spring forage indices

We can compare the indices relevant to bluefish on the same scale: inshore (Albatross strata), inshore bluefish, offshore bluefish, and further out, plus the state and federal waters split (Fig. 3.9).

Figure 3.9: Bluefish-related Spring forage indices

Proportions of prey in each bluefish relevant area for spring can be compared (here as a proportion of the full MABGB index; Fig. 3.10).

Figure 3.10: Spring forage indices as proportion of the Mid Atlantic and Georges Bank area.

3.3.5 Spring predicted ln-density

The VAST model predicts density of forage fish across the entire model domain for each year (Fig. 3.11).

Figure 3.11: Yearly maps of VAST model estimated forage fish density for Spring (months 1-6).

3.3.6 Spring Diagnostics

Diagnostics shown for the spring model are as described above for the fall model.

3.3.6.1 Data_and_knots

Spring diagnostics:Data_and_knots

3.3.6.2 Data_by_year

Spring diagnostics:Data_by_year

3.3.6.3 quantile_residuals

Spring diagnostics:quantile_residuals

3.3.6.4 quantile_residuals_on_map

Spring diagnostics:quantile_residuals_on_map

3.3.6.5 Aniso

Spring diagnostics:Aniso

3.3.6.6 center_of_gravity

Spring diagnostics:center_of_gravity

3.3.7 Annual Index

The Annual forage index uses data from all months each year (1-12). The same plots are presented as those above for Fall and Spring indices. Fig. 3.12 shows all annual forage index time series with standard errors.

Figure 3.12: All Annual forage indices

We can compare the annual indices relevant to bluefish on the same scale: inshore (Albatross strata), inshore bluefish, offshore bluefish, and further out, plus the state and federal waters split (Fig. 3.13).

Figure 3.13: Bluefish-related Annual forage indices

Proportions of prey in each bluefish relevant area for the annual model can be compared (here as a proportion of the full MABGB index; Fig. 3.14).

Figure 3.14: Annual forage indices as proportion of the Mid Atlantic and Georges Bank area.

3.3.8 Annual predicted ln-density

The VAST model predicts density of forage fish across the entire model domain for each year (Fig. 3.15).

Figure 3.15: Yearly maps of VAST model estimated annual forage fish density (months 1-12).

3.3.9 Annual Diagnostics

Diagnostics shown for the annual model are as described above for the falland spring models.

3.3.9.1 Data_and_knots

Annual diagnostics:Data_and_knots

3.3.9.2 Data_by_year

Annual diagnostics:Data_by_year

3.3.9.3 quantile_residuals

Annual diagnostics:quantile_residuals

3.3.9.4 quantile_residuals_on_map

Annual diagnostics:quantile_residuals_on_map

3.3.9.5 Aniso

Annual diagnostics:Aniso

3.3.9.6 center_of_gravity

Annual diagnostics:center_of_gravity

The full results of all models are archived on google drive rather than github to save space.

3.4 Woods Hole stock assessment model (WHAM) input time series

Current inputs for the bluefish assessment implemented in WHAM (Stock and Miller 2021) include a subset of the stratum-specific indices calculated above. The index is calculated for each season and the annual model for several regions relevant to the bluefish assessment:

1. Albatross New (AlbNew) includes all inshore and new offshore survey strata (largest area, combines area in Figs. 2.3, 2.4, and 2.5)
2. Albatross Old (AlbOld) includes all inshore survey strata (combines area in Figs. 2.3 and 2.4)
3. Bigelow New (BigNew) includes the subset of inshore survey strata that can be sampled by the R/V Henry Bigelow plus new offshore strata (combines area in Figs. 2.4 and 2.5)
4. Bigelow Old (BigOld) includes the subset of inshore survey strata that can be sampled by the R/V Henry Bigelow (area in Fig. 2.4)
5. StateWaters includes the coastline to 3 nautical miles offshore (smallest area, in Fig. 2.6)

WHAM model inputs were processed with the WHAMinputs.R script.

Indices in the strata above for input in to WHAM are plotted below (Figs. 3.16, 3.17, 3.18).

3.4.1 WHAM forage index time series Fall

Figure 3.16: Time series of VAST estimated fall forage indices for input into WHAM, 1985-2021

3.4.2 WHAM forage index time series Spring

Figure 3.17: Time series of VAST estimated spring forage indices for input into WHAM, 1985-2021

3.4.3 WHAM forage index time series Annual

Figure 3.18: Time series of VAST estimated annual forage indices for input into WHAM, 1985-2021

All of the above indices were provided for testing within the WHAM-based bluefish assessment model.

4 Discussion

Bluefish are generalist predators. Understanding how bluefish availability may be affected by their prey field means looking across many potential prey species, including species not well sampled by bottom trawl surveys. We characterized aggregate forage fish trends using piscivore stomach data as input observations to a vector autoregressive spatio-temporal (VAST) model. The resulting model-derived forage indices suggest that aggregate forage fish biomass has fluctuated in both space and time in the Northeast US, especially in areas relevant to the bluefish assessment. Evaluating these indices within the bluefish assessment model as covariates reflecting bluefish availability to fisheries and surveys is in progress.

Although time series of aggregate forage biomass for direct comparison with this work are lacking, we think, due to good sample sizes and VAST model diagnostics, that this is the best representation we have of aggregate forage fish trends in this region. Fish stomach data has been used to index prey abundance, especially for data-poor prey taxa poorly sampled by bottom trawls, in multiple ecosystems worldwide (Fahrig et al. 1993; Link 2004; Mills et al. 2007; Cook and Bundy 2012; Lasley-Rasher et al. 2015; Rohan and Buckley 2017; Sydeman et al. 2022). Integrating fish stomach content information within a VAST model to index prey biomass was shown to be successful with selected predators of Atlantic herring in the Northeast US (Ng et al. 2021), using a subset of the data we used here. In particular, reasonable agreement was found between the Atlantic cod diet-derived Atlantic herring biomass index and Atlantic herring assessment biomass trends. However, herring indices based on stomach contents from other individual predators (e.g., silver hake) had poor agreement with herring assessment trends (Ng et al. 2021). Ng et al. (2021) noted that decadal trends generally agreed across predators between diet-based and stock assessment Atlantic herring biomass indices, but that shorter term trends varied considerably by seasons and predators. Further, changes in spatial overlap between individual predators and prey over time was noted by Ng et al. (2021) as a potential issue with using diet-based indices of Atlantic herring abundance.

Aggregation of predators and seasons was one recommendation for future work from Ng et al. (2021) to address some of the issues with using predator stomach data to evaluate prey biomass. Our approach combining stomach data across multiple piscivore predators was intended to improve sample size for the aggregate prey index so that our results integrate across individual predator sampling variability and changes in predator-prey overlap, and perhaps clarify the longer term, decadal trends in forage fish biomass. Further, evaluating forage species in aggregate should also reduce issues of changing overlap of predators with individual prey species. However, we don’t have an aggregate forage fish assessment time series for comparison. Using survey biomass catch of forage fish directly can be problematic, including issues with poor sample size or missing data for key small pelagic fish (anchovies, sandlance) and cephalopods (Mills et al. 2007; Rohan and Buckley 2017). Further, bottom trawl survey biomass time series for assessed forage species can be at odds with other assessment inputs and assessment-estimated biomass trends in the Northeast US. For example, Atlantic mackerel bottom trawl survey indices have generally increased while assessment estimated biomass, largely driven by an egg survey and information on catch, is decreasing (NEFSC 2018). Therefore, our forage index represents an alternative, possibly more complete assessment of aggregate forage fish in this region than can be obtained directly from bottom trawl surveys.

Changes in forage fish are of interest across seasons, but changes in the fall index nearshore may be most closely aligned with abundance indices used in the bluefish assessment. Bluefish migrate into northern coastal waters in spring and return south in late fall to overwinter (Collette and Klein-Macphee 2002). Although bluefish catch patterns vary by state, recreational fishing dominates coastwide bluefish landings, and most recreational fishing activity (nearly 70% of landings) takes place in July-October (2022 bluefish research track summary report), months included in the fall forage index (months 7-12). Further, northern states, included in the spatial extent of the forage index, see peak recreational fishing during these months. Therefore, our fall forage indices may provide a reasonable match to the recreational catch per unit effort index used in the assessment. Similarly, fall forage indices for fishery independent survey strata also temporally and spatially align with both the fall NEFSC bottom trawl survey and fall NEAMAP bottom trawl survey stratified mean catch-per-tow used in the bluefish assessment (2022 bluefish research track summary report).

While the fall forage fish indices are temporally aligned with bluefish assessment inputs and spatially aligned with two trawl survey indices used in the assessment, improvements in spatial overlap with recreational fisheries and other survey indices could be considered in the future. A large proportion (51%) of bluefish recreational landings come from inland waters: bays and estuaries including Chesapeake Bay, Delaware Bay, and Long Island Sound, with the next largest proportion (42%) coming from state waters extending 3 nautical miles from the coastline (2022 bluefish research track summary report). The current forage index does not cover inland waters, aside from Narragansett Bay and Buzzards Bay. Diet data are available for Chesapeake Bay from the ChesMMAP survey, which could be added to the VAST model in the future. Less diet information is available for the portion of the bluefish range south of Cape Hatteras, although some collections have taken place. Investigation of sources of diet information, or possibly direct forage fish surveys for inland and southern areas would be worthwhile to see whether data are adequate to cover the full range of bluefish. In addition, it is worth exploring whether including indices as area-expanded estimates or as unexpanded forage density in each area of interest provides more information to the bluefish assessment model.

The general method of using stomach information to estimate prey indices can be extended to address other predators with different prey fields (e.g., benthic feeders; Link (2004); Lasley-Rasher et al. (2015); Rohan and Buckley (2017)) or different spatial distributions. Further adjustments to the VAST model could explore the sensitivity of aggregate forage trends to the inclusion or exclusion of predators from diet sampling over time and the level of species identification in diets. For example, we aggregated predators most similar to bluefish based on analysis of the full diet database to minimize impacts from uneven sampling of individual predators, but some differences in predator sampling over time (e.g. squid stomachs sampled up to 2000s, not thereafter) may still affect the forage index. More importantly, we had to leave out the “fish unidentified” prey category which comprised 25-30% of bluefish diet because while these prey may be forage fish for bluefish, there is no guarantee they would be for other sampled predators we included in the model. Increased use of a broad “fish unidentified” category would degrade the use of this index in the future, so continued investment in high quality stomach data collection is crucial. In addition, including more information on the process of fish predation (which we somewhat controlled for as catchability covariates for number of predators, predator size, and temperature) may help to refine aggregate forage indices. For example, in this initial model we have not accounted for functional responses of predators, but new research may allow us to do so in future iterations (Smith and Smith 2020; Robertson et al. 2022).

Overall, this index provides insight into temporal and spatial variation at the forage fish community level, which is important both for individual predators and for ecosystem assessment. Forage fish link lower trophic level productivity with larger fish important for human consumption and recreation as well as for protected species, so understanding aggregate forage dynamics within an ecosystem may support analysis related to management for multiple living resources (Smith et al. 2011; Essington et al. 2015; Levin et al. 2016; Punt et al. 2016; Tommasi et al. 2017; Soudijn et al. 2021). While aggregate indices of forage fish may be inherently more stable than individual forage species population trajectories, we still find substantial fluctuations in this forage index. Investigation into drivers of forage fish (and predator) spatial and temporal shifts demonstrate substantial variability. Many factors influence forage distribution and abundance, including environmental drivers changing habitat and impacting species differently, resulting in often unclear or mixed signals across taxa, although general trends in the Northeast US are towards the northeast and into deeper water (Fredston et al. 2021; Suca et al. 2021). This initial aggregate forage index provides the opportunity to investigate whether temporal and spatial trends are coherent with aggregate zooplankton indices and or spatial and temporal patterns in environmental drivers. Ongoing analyses of ecosystem linkages with the forage index may provide insight both for improving the index for future predator stock assessments and for ecosystem reporting in the Northeast US (NEFSC 2022).

References

Collette, B. B., and G. Klein-Macphee. 2002. Bigelow and Schroeder’s Fishes of the Gulf of Maine, Third Edition3rd ed. edition. Smithsonian Books, Washington, DC.

Cook, A. M., and A. Bundy. 2012. Use of fishes as sampling tools for understanding biodiversity and ecosystem functioning in the ocean. Marine Ecology Progress Series 454:1–18.

Essington, T. E., P. E. Moriarty, H. E. Froehlich, E. E. Hodgson, L. E. Koehn, K. L. Oken, M. C. Siple, and C. C. Stawitz. 2015. Fishing amplifies forage fish population collapses. Proceedings of the National Academy of Sciences 112(21):6648–6652.

Fahrig, L., G. R. Lilly, and D. S. Miller. 1993. Predator Stomachs as Sampling Tools for Prey Distribution: Atlantic Cod ( Gadus morhua ) and Capelin ( Mallotus villosus ). Canadian Journal of Fisheries and Aquatic Sciences 50(7):1541–1547.

Fredston, A., M. Pinsky, R. L. Selden, C. Szuwalski, J. T. Thorson, S. D. Gaines, and B. S. Halpern. 2021. Range edges of North American marine species are tracking temperature over decades. Global Change Biology 27(13):3145–3156.

Garrison, L., and J. Link. 2000. Dietary guild structure of the fish community in the Northeast United States continental shelf ecosystem. Marine Ecology Progress Series 202:231–240.

Lasley-Rasher, R. S., D. C. Brady, B. E. Smith, and P. A. Jumars. 2015. It takes guts to locate elusive crustacean prey. Marine Ecology Progress Series 538:1–12.

Levin, P. S., T. B. Francis, and N. G. Taylor. 2016. Thirty-two essential questions for understanding the social–ecological system of forage fish: The case of pacific herring. Ecosystem Health and Sustainability 2(4):e01213.

Link, J. S. 2004. Using fish stomachs as samplers of the benthos: Integrating long-term and broad scales. Marine Ecology Progress Series 269:265–275.

Mills, K. L., T. Laidig, S. Ralston, and W. J. Sydeman. 2007. Diets of top predators indicate pelagic juvenile rockfish (Sebastes spp.) Abundance in the California Current System. Fisheries Oceanography 16(3):273–283.

NEFSC. 2018. 64th Northeast Regional Stock Assessment Workshop (64th SAW) Assessment Report. Northeast Fisheries Science Center (U.S.) reference document.

NEFSC. 2022. 2022 State of the Ecosystem Mid-Atlantic. Northeast Fisheries Science Center (U.S.):48.

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

Northeast Area Monitoring & Assessment Program (NEAMAP), V. S. M. & Assessment Program (VASMAP), R. J. Latour, J. Gartland, and C. F. Bonzek. 2021. Monitoring Living Marine Resources in the Mid‐Atlantic Bight: Virginia Institute of Marine Science, William & Mary.

Perretti, C. T., and J. T. Thorson. 2019. Spatio-temporal dynamics of summer flounder (Paralichthys dentatus) on the Northeast US shelf. Fisheries Research 215:62–68.

Punt, A. E., A. D. MacCall, T. E. Essington, T. B. Francis, F. Hurtado-Ferro, K. F. Johnson, I. C. Kaplan, L. E. Koehn, P. S. Levin, and W. J. Sydeman. 2016. Exploring the implications of the harvest control rule for Pacific sardine, accounting for predator dynamics: A MICE model. Ecological Modelling 337:79–95.

Reynolds, R. W., T. M. Smith, C. Liu, D. B. Chelton, K. S. Casey, and M. G. Schlax. 2007. Daily High-Resolution-Blended Analyses for Sea Surface Temperature. Journal of Climate 20(22):5473–5496.

Robertson, M. D., M. Koen-Alonso, P. M. Regular, N. Cadigan, and F. Zhang. 2022. Accounting for a nonlinear functional response when estimating prey dynamics using predator diet data. Methods in Ecology and Evolution 13(4):880–893.

Rohan, S. K., and T. W. Buckley. 2017. Spatial and ontogenetic patterns of Pacific cod (Gadus macrocephalus Tilesius) predation on octopus in the eastern Bering Sea. Environmental Biology of Fishes 100(4):361–373.

Saha, S., Korak, X. Zhao, H. Zhang, K. S. Casey, D. Zhang, S. Baker-Yeboah, K. A. Kilpatrick, R. H. Evans, T. Ryan, and J. M. Relph. 2018. AVHRR Pathfinder version 5.3 level 3 collated (L3C) global 4km sea surface temperature for 1981-Present. [1985-2021]. NOAA National Centers for Environmental Information. Dataset.

Schoener, T. W. 1970. Nonsynchronous Spatial Overlap of Lizards in Patchy Habitats. Ecology 51(3):408–418.

Smith, A. D., C. J. Brown, C. M. Bulman, E. A. Fulton, P. Johnson, I. C. Kaplan, H. Lozano-Montes, S. Mackinson, M. Marzloff, L. J. Shannon, Y. J. Shin, and J. Tam. 2011. Impacts of fishing low-trophic level species on marine ecosystems. Science 333(6046):1147–50.

Smith, B. E., and J. S. Link. 2010. The Trophic Dynamics of 50 Finfish and 2 Squid Species on the Northeast US Continental Shelf. NOAA Technichal Memorandum NMFS-NE-216. National Marine Fisheries Service, 166 Water Street, Woods Hole, MA 02543-1026.

Smith, B. E., and L. A. Smith. 2020. Multispecies functional responses reveal reduced predation at high prey densities and varied responses among and within trophic groups. Fish and Fisheries 21(5):891–905.

Soudijn, F. H., P. Daniël van Denderen, M. Heino, U. Dieckmann, and A. M. de Roos. 2021. Harvesting forage fish can prevent fishing-induced population collapses of large piscivorous fish. Proceedings of the National Academy of Sciences 118(6):e1917079118.

Stock, B. C., and T. J. Miller. 2021. The Woods Hole Assessment Model (WHAM): A general state-space assessment framework that incorporates time- and age-varying processes via random effects and links to environmental covariates. Fisheries Research 240:105967.

Suca, J. J., J. J. Deroba, D. E. Richardson, R. Ji, and J. K. Llopiz. 2021. Environmental drivers and trends in forage fish occupancy of the Northeast US shelf. ICES Journal of Marine Science (fsab214).

Sydeman, W. J., S. A. Thompson, J. F. Piatt, S. G. Zador, and M. W. Dorn. 2022. Integrating seabird dietary and groundfish stock assessment data: Can puffins predict pollock spawning stock biomass in the North Pacific? Fish and Fisheries 23(1):213–226.

Team, R. C. 2021. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.

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

Thorson, J. T., and L. A. K. Barnett. 2017. Comparing estimates of abundance trends and distribution shifts using single- and multispecies models of fishes and biogenic habitat. ICES Journal of Marine Science 74(5):1311–1321.

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

Tommasi, D., C. A. Stock, K. Pegion, G. A. Vecchi, R. D. Methot, M. A. Alexander, and D. M. Checkley. 2017. Improved management of small pelagic fisheries through seasonal climate prediction. Ecological Applications 27(2):378–388.

Zuur, A. F., E. N. Ieno, N. Walker, A. A. Saveliev, and G. M. Smith. 2009. Mixed effects models and extensions in ecology with R. Springer, New York, NY.