Bluefish forage index thoughts
Sarah Gaichas
10 November, 2021
Does prey drive availability of bluefish?
Proposal: create a “prey index” to evaluate changes in bluefish prey over time and in space using VAST
Pattern on (Ng et al., 2021), which used predator stomach data to create a biomass index for a single prey, herring. Expected biomass of herring per stomach was estimated with 2 linear predictors, n herring per stomach and av wt of herring in a stomach.
However, we want a biomass index for “bluefish prey” rather than a single prey species. Further, we want to include inshore and offshore regions by combining surveys as was done for summer flounder biomass in (Perretti and Thorson, 2019). Finally, since bluefish themselves are somewhat sparsely sampled by the surveys, I think we should consider aggregating all predators that have a similar diet composition to bluefish to better represent bluefish prey biomass.
I think what we want is to characterize weight of bluefish prey from all piscivores caught at each location and model that over time/space.
Covariates explaining patterns in bluefish prey index could include number of predators, species composition of predators, size composition of predators at each location? and could throw in temperature.
What are 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 here for NEFSC survey and here for NEAMAP survey bluefish diet composition summaries.
# use same prey list as from shiny summaries
bluefishaggdiet <- read_csv(here("datfromshiny","DC.Bluefish (Pomatomus saltatrix).All.2021-10-07.csv"))
# decadal diet
diet70 <- read_csv(here("datfromshiny","DC.Bluefish (Pomatomus saltatrix).1970s.2021-10-07.csv")) %>%
mutate(Decade = 1970)
diet80 <- read_csv(here("datfromshiny","DC.Bluefish (Pomatomus saltatrix).1980s.2021-10-07.csv"))%>%
mutate(Decade = 1980)
diet90 <- read_csv(here("datfromshiny","DC.Bluefish (Pomatomus saltatrix).1990s.2021-10-07.csv"))%>%
mutate(Decade = 1990)
diet00 <- read_csv(here("datfromshiny","DC.Bluefish (Pomatomus saltatrix).2000s.2021-10-07.csv"))%>%
mutate(Decade = 2000)
diet10 <- read_csv(here("datfromshiny","DC.Bluefish (Pomatomus saltatrix).2010s.2021-10-07.csv"))%>%
mutate(Decade = 2010)
bluefishdecadediet <- bind_rows(diet70, diet80, diet90, diet00, diet10)
# seasonal diet
dietspring <- read_csv(here("datfromshiny","DC.Bluefish (Pomatomus saltatrix).SPRING.2021-10-07.csv"))%>%
mutate(Season = "Spring")
dietfall <- read_csv(here("datfromshiny","DC.Bluefish (Pomatomus saltatrix).FALL.2021-10-07.csv"))%>%
mutate(Season = "Fall")
bluefishseasondiet <- bind_rows(dietspring, dietfall)
# regional diet
dietMAB <- read_csv(here("datfromshiny","DC.Bluefish (Pomatomus saltatrix).MID-ATLANTIC BIGHT.2021-10-07.csv")) %>%
mutate(Region = "MAB")
dietSNE <- read_csv(here("datfromshiny","DC.Bluefish (Pomatomus saltatrix).SOUTHERN NEW ENGLAND.2021-10-07.csv"))%>%
mutate(Region = "SNE")
dietGB <- read_csv(here("datfromshiny","DC.Bluefish (Pomatomus saltatrix).GEORGES BANK.2021-10-07.csv"))%>%
mutate(Region = "GB")
bluefishregiondiet <- bind_rows(dietMAB, dietSNE, dietGB)
# bluefish size diet
dietSM <- read_csv(here("datfromshiny","DC.Bluefish (Pomatomus saltatrix).S.2021-10-07.csv")) %>%
mutate(Size = "Small")
dietMED <- read_csv(here("datfromshiny","DC.Bluefish (Pomatomus saltatrix).M.2021-10-07.csv"))%>%
mutate(Size = "Medium")
dietLG <- read_csv(here("datfromshiny","DC.Bluefish (Pomatomus saltatrix).L.2021-10-07.csv"))%>%
mutate(Size = "Large")
bluefishsizediet <- bind_rows(dietSM, dietMED, dietLG)
preylist <- unique(c(bluefishaggdiet$Prey, bluefishdecadediet$Prey, bluefishregiondiet$Prey, bluefishseasondiet$Prey, bluefishsizediet$Prey))35 species from preylist:
Unidentified fish, Engraulidae, Anchoa mitchilli, Peprilus triacanthus, Loligo sp, Ammodytes sp, Pomatomus saltatrix, Clupea harengus, Merluccius bilinearis, Anchoa hepsetus, Cephalopoda, Etrumeus teres, Clupeidae, Illex sp, Stenotomus chrysops, Brevoortia tyrannus, Alosa pseudoharengus, Macrozoarces americanus, Loligo pealeii, Pleuronectiformes, Prionotus sp, Scomberesox saurus, Urophycis chuss, Alosa aestivalis, Melanogrammus aeglefinus, Scomber sp, Ctenophora, Engraulis eurystole, Scomber scombrus, WDP, Myoxocephalus octodecemspinosus, Euphausiacea, Urophycis regia, Squalus acanthias, Scophthalmus aquosus
We can aggregate further, but are there any in particular we want to track? Are any categories more important?
Has the total amount of bluefish prey changed over time?
Because bluefish prey are mainly small pelagic species not well sampled by trawl surveys, direct estimates of prey biomass from trawl surveys are unlikely to be useful to understand changes in prey biomass over time. However, we could still try this in aggregate if we have time.
Multiple predators that are more effectively sampled by trawl surveys also consume bluefish prey. Therefore, using a suite of predators as “samplers” of bluefish prey may more effectively estimate changes in prey over time. This method has been applied to develop an index of Atlantic herring in the Northeast US using the NEFSC bottom trawl survey diet data (Ng et al., 2021). In our case, we will want to apply the method to multiple prey rather than a single prey species.
It may be useful to track prey categories from predator diets (family level aggregation) rather than species level or fully aggregated “prey”. We can evaluate how much data is available at each level of aggregation.
Hypothesis: we would not expect significant changes over time in an index aggregate prey biomass over bluefish’s range; they prey on a wide variety of species. However, we may see fluctuations in prey types that could be informative.
Caveat: the surveys here do not cover bluefish’s range (sparse south of Cape Hatteras).
Has the spatial distribution of bluefish prey changed over time?
To evaluate whether bluefish prey distribution has changed, we should (at a minimum) combine nearshore and offshore surveys in as many seasons as feasible. The NEFSC and NEAMAP bottom trawl surveys both operate in spring and fall, and have been combined to analyze changes in biomass and distribution of summer flounder (Perretti and Thorson, 2019).
Methods
We will use VAST (Thorson and Barnett, 2017; Thorson, 2019) to evaluate changes in bluefish prey biomass and distribution over time.
\[\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) \\ + \nu_1(c_i, t_i) + \zeta_1(i) - \iota(c_i, t_i)\]
Decisions (15 major)
Spatial domain
Bluefish range is coastwide but these surveys cover only northern half. Could limit domain to Hatteras north to Canada on continental shelf, span of NEFSC and NEAMAP bottom trawl surveys. Therefore we would have a “northeast US bluefish prey index.”
Categories (Species/sizes)
Which predators to include?
Fish feeding guilds based on dietary overlap have been defined for the Northeast US shelf based on NEFSC trawl survey diet data from 1973-1997 (Garrison and Link, 2000). In this analysis, all size classes of bluefish were classified as piscivores. A reasonable first cut would be to use all predator/size combinations identifed as piscivores in Garrison and Link (2000). This would include:
# code for piscivore table
piscivores <- ecodata::species_groupings %>%
select(COMNAME, SizeCat, Garrison.Link) %>%
filter(!is.na(Garrison.Link),
Garrison.Link == "Piscivores") %>%
mutate(PiscGuild = case_when(COMNAME == "WINTER SKATE" ~ "c",
COMNAME == "WEAKFISH" ~ "b",
COMNAME == "BLUEFISH" & SizeCat == "S" ~ "b",
TRUE ~ "a")) %>%
distinct()
datatable(piscivores, rownames = FALSE,
caption = "Piscivores feeding guild from Garrison and Link 2000",
options = list(pageLength = 25))Should predators be aggregated prior to analysis or kept separate?
Options are each predator and size class in the piscivore guild separate, each predator separate (aggregate size), each piscivore guild separate (“a”, “b”, and “c”), or all piscivores combined. A model with only bluefish predators could also be done.
Should we account for predator size in an aggregation? Garrison and Link (2000) noted that small bluefish along with weakfish formed a subset of piscivores with predominantly anchovy prey (PiscGuild “b” in the table). Large summer flounder and large sharpnose sharks had the highest similarity to medium and large bluefish, but were not distinguished as a subset of the guild from the remaining piscivores (Garrison and Link, 2000). Winter skates were a separate piscivore guild “c” so we could keep them or not.
There have been 20+ more years of data collected since this analysis. It seems unlikely that this list would change much, but if we had Lance’s code we could possibly rerun? Maybe Brian has?
Is it useful for the assessment to characterize prey of small (up to 30 cm) bluefish separately from 31+ cm bluefish? If so we could separate these predator guilds. There may be more weakfish diet data in the NEAMAP survey than in the NEFSC survey, so joining the surveys will be important if we want to do this.
Combine predators by ignoring species or treat each as a “survey” contributing to prey index estimate? Do we need pred biomass (more like predator-expanded-stomach-contents, (Grüss et al., 2020)) if we treat them as a survey? I don’t think so as this is getting at predator diet composition and absolute consumption, while we are interested in prey biomass itself.
Which prey to include?
Entire list of bluefish prey? 35 categories from NEFSC summaries. In the test dataset there are over a hundred individual prey categories.
How to aggregate prey?
Taxonomic? Size? Habitat? Combination?
Avoid modeling all 35+ identified prey separately! Taxonomic aggregation works well with how species id is done in stomach samples as many are family level or above.
Bay anchovy separate? All other anchovies, herrings/menhaden separate? squids
Which variable –> Which model structure
Presence/absence of prey, count of prey, biomass? Prey count is not practical as partial remains are common. Prey biomass based on stomach contents. Delta model. Poisson-link model alternative which has been used in both diet examples.
Univariate: all prey together, all pisivores together with piscivore attributes (n spp, mean length, spp diversity, one key pred dominant, ???) as factor(s)? Density covariates? Catchability covariates?
Univariate complex: each prey category modeled, all piscivores together
Multivariate: all prey categories modeled together, all piscivores together
(Could try a model with only bluefish as a predator but unsure it would work)
Seasons (Spring, Fall) separate for each analysis.
VAST questions
Explain factors (relates multivariate responses) vs density covariates (explains spatial and temporal patterns due to habitat) vs catchability covariates (predator attributes seem to fit this for sampling prey, but do we need to account for survey catchability of predators? due to observation process). Can any of these be the same thing?
Advice for specifying/separating the above
We are looking maybe at a hierarchical structure: Survey catches predators (model accounts for observation process interacting with ecology of predators) Predators sample prey (model accounts for interacting ecology of predators and prey) We want to reconstruct the prey
Data
NEFSC bottom trawl diet data
Initial test food habits dataset from 2018 ECSA has food habits 1973 to 2016. Will get fuller up to date info from Brian Smith shortly.
# code for joining food habits data to piscivore table
load(here("fhdat", "allfhsg_2016.RData"))
fh.nefsc.pisc <- allfhsg %>%
filter(pynam != "EMPTY") %>%
left_join(piscivores, by = c("pdcomnam" = "COMNAME",
"sizecat" = "SizeCat")) %>%
filter(!is.na(Garrison.Link))Here are size distributions of piscivores in my trial NEFSC diet dataset 1973 to 2016. We don’t have weakfish in this dataset anyway:
# code for making figure
fh.nefsc.pisc %>%
group_by(year, season, station, pdid) %>%
ggplot(aes(x=pdlen, fill=pdcomnam)) +
geom_histogram(stat = "count") +
theme_minimal() +
theme(legend.position = "none") +
facet_wrap(~pdcomnam, scales = "free")
Length distribution of piscivore predators in dataset.
We can get down to 17 prey categories (“analsci”) for the piscivore dataset by filtering using the most common prey identified in the summaries:
fh.nefsc.pisc.bfprey <- fh.nefsc.pisc %>%
filter(tolower(pynam) %in% tolower(preylist))17 categories using analsci column
CEPHALOPODA, COTTIDAE, PLEURONECTIFORMES, CLUPEIDAE, GADIDAE, ZOARCIDAE, STROMATEIDAE, POMATOMIDAE, SCOMBRIDAE, BOTHIDAE, CTENOPHORA, AMMODYTIDAE, SPARIDAE, ENGRAULIDAE, SHARK, TRIGLIDAE, OTHER FISH
So instead of averaging prey weights in stomachs per species at a station, I should average all over all piscivore stomachs at a station. Then do an index for each prey category from the aggreate piscivores.
So add a column mapping pynam to blueprey and filter by that. Then sum that preywt at each station and divide by n preds using the piscivore dataset. Then that is the input? average prey biomass per stomach per station to get number per stomach and average weight per stomach, ultimately estimating expected biomass of prey per piscivore stomach?
And what are the covariates? Could be number of predators, species composition of predators, size composition of predators at each location? and could throw in temperature.