Update Input Data to remove modeled prey

The bluefish forage index includes the most common prey of bluefish, including managed forage species such as Atlantic mackerel and Atlantic herring.

Our challenge is to create a forage index for the combined predators modeled in Hydra, that does not include the combined prey modeled in Hydra.

See the mskeyrun overview at https://noaa-edab.github.io/ms-keyrun/articles/mskeyrun.html for the full list.

Hydra modeled predators: Atlantic cod (Gadus morhua), Goosefish (Lophius americanus), Silver hake (Merluccius bilinearis), Spiny dogfish (Squalus acanthias)

Hydra modeled prey: Atlantic herring (Clupea harengus), Atlantic mackerel (Scomber scombrus), (and Silver hake)

The Hydra modeled predators all clustered with bluefish by diet similarity. Therefore, I think we can keep the same list of piscivores to sample the prey field.

We don’t want the non-pelagic prey of these predators in the index, only other pelagic forage.

So we can resort the prey list of the piscivores to see if we would add any other prey and then remove Atlantic mackerel and Atlantic herring from the prey list.

We checked the NEAMAP prey list and found we didn’t include silver hake and Atlantic mackerel due to the 10 stomach cutoff, and Atlantic herring are mixed with Clupeids unident which include higher proportions of river herrings that we want to keep. Therefore we have made no changes to the NEAMAP data at present.

[However, we may change both inputs in the future to include all species that made the 10 stomach cutoff in either survey. Stay tuned. It doesn’t change the input for this purpose though.]

Working group decision after the February 2022 Data Meeting was to update the predator list based on the most recent diet similarity matrix by Smith posted on the NEFSC shiny app, which results in a new Piscivore guild relative to our initial models based on Garrison and Link (2000).

We investigated alternative cluster algorithms here.

Input NEFSC food habits overlap matrix:

dietoverlap <- read_csv(here("datfromshiny/tgmat.2022-02-15.csv"))

Get NEFSC food habits data (1973-2020 combined with 2021):

# Oct 21 2022: save it for posterity, use in case Laurel's repo updates
#save(allfh, file = "fhdat/allfh.Rdata")
load(here("fhdat/allfh.Rdata"))

# October 8 2022: add NEFSC 2021 data
load(here("fhdat/allfh21.Rdata"))

allfh <- allfh %>%
  dplyr::bind_rows(allfh21)

Generate the Piscivore list based on clustering with the “complete” algorithm. Ensure our Hydra predators are in the cluster.

# follows example here https://cran.r-project.org/web/packages/dendextend/vignettes/Cluster_Analysis.html

library(dendextend)

d_dietoverlap <- dist(dietoverlap)

guilds <- hclust(d_dietoverlap)

#plot(guilds)

dend <- as.dendrogram(guilds)

dend <- dendextend::rotate(dend, 1:136)

dend <- color_branches(dend, k=6) # Brian uses 6 categories

labels(dend) <- paste(as.character(names(dietoverlap[-1]))[order.dendrogram(dend)],
                           "(",labels(dend),")", 
                           sep = "")

dend <- hang.dendrogram(dend,hang_height=0.1)

# reduce the size of the labels:
# dend <- assign_values_to_leaves_nodePar(dend, 0.5, "lab.cex")
dend <- set(dend, "labels_cex", 0.75)
# And plot:
# par(mar = c(3,3,3,7))
# plot(dend, 
#      main = "Clustered NEFSC diet data, (complete)
#      (the labels give the predator species/size)", 
#      horiz =  TRUE,  nodePar = list(cex = .007))
#legend("topleft", legend = iris_species, fill = rainbow_hcl(3))

circlize::circos.par(start.degree = 90)
par(mar = rep(0,4))
circlize_dendrogram(dend,
                    labels_track_height=0.3,
                    dend_track_height = 0.6)

This is the list of predators and sizes:

# list of species in node with all three bluefish sizes
pisccomplete <- partition_leaves(dend)[[
  which_node(dend, c("Bluefish..S(37)", "Bluefish..M(36)", "Bluefish..L(35)"))
]]

pisccomplete
##  [1] "Striped bass..M(106)"           "Sea raven..L(89)"              
##  [3] "Striped bass..L(105)"           "Bluefish..S(37)"               
##  [5] "Weakfish..M(117)"               "Fourspot flounder..L(50)"      
##  [7] "Northern shortfin squid..M(71)" "Longfin squid..M(63)"          
##  [9] "Longfin squid..S(64)"           "Northern shortfin squid..S(72)"
## [11] "Pollock..L(78)"                 "White hake..M(120)"            
## [13] "Red hake..L(82)"                "Silver hake..M(93)"            
## [15] "Atlantic halibut..M(17)"        "Pollock..XL(81)"               
## [17] "White hake..L(119)"             "Silver hake..L(92)"            
## [19] "Cusk..L(45)"                    "Goosefish..XL(56)"             
## [21] "Bluefish..L(35)"                "Goosefish..L(53)"              
## [23] "Goosefish..M(54)"               "Goosefish..S(55)"              
## [25] "Spiny dogfish..M(100)"          "Thorny skate..XL(116)"         
## [27] "Atlantic halibut..L(16)"        "Sea raven..M(90)"              
## [29] "Sea raven..S(91)"               "Atlantic cod..XL(13)"          
## [31] "Spiny dogfish..L(99)"           "Spotted hake..M(103)"          
## [33] "Summer flounder..M(111)"        "Summer flounder..L(110)"       
## [35] "Bluefish..M(36)"                "Buckler dory..M(38)"

We have all sizes of goosefish (monkfish), XL Atlantic cod, M and L Spiny dogfish, M and L Silver hake.

Lets call this close enough to define other forage.

Create the NEFSC input dataset.

First filter to the predators and identify the bluefish prey:

pisccompletedf <- data.frame("COMNAME" = toupper(str_remove(pisccomplete, "\\..*")),
                              "SizeCat" = str_remove(str_extract(pisccomplete, "\\..*[:upper:]+"), "\\.."),
                              "feedguild" = "pisccomplete")
 
 fh.nefsc.pisc.pisccomplete <- allfh %>%
  #filter(pynam != "EMPTY") %>%
  left_join(pisccompletedf, by = c("pdcomnam" = "COMNAME",
                               "sizecat" = "SizeCat")) %>%
  filter(!is.na(feedguild))
 
preycount  <- fh.nefsc.pisc.pisccomplete %>%
   #group_by(year, season, pdcomnam, pynam) %>%
   group_by(pdcomnam, pynam) %>%
   summarise(count = n()) %>%
   #arrange(desc(count))
   pivot_wider(names_from = pdcomnam, values_from = count) %>%
   select(pynam, "ATLANTIC COD", "BLUEFISH", "GOOSEFISH", "SILVER HAKE", "SPINY DOGFISH")

datatable(preycount)

Remove modeled prey

Now removing the three modeled species from the prey list to create new data files with “_unmod” appended. This code removes “MERLUCCIUS BILINEARIS”, “SCOMBER SCOMBRUS”, “CLUPEA HARENGUS”, and “CLUPEIDAE” from the NEFSC dataset, as well as the useless categories empty, blown, fish, animal remains, etc., and the large family level invertebrate categories.

gencomlist <- allfh %>%
  select(pynam, gencom2) %>%
  distinct()

blueprey.unmod <- preycount %>% 
  filter(BLUEFISH > 9) %>%
  filter(!pynam %in% c("EMPTY", "BLOWN",
                       "FISH", "OSTEICHTHYES",
                       "ANIMAL REMAINS",
                       "FISH SCALES",
                       "MERLUCCIUS BILINEARIS",  # remove modeled species
                       "SCOMBER SCOMBRUS",  # remove modeled species
                       "CLUPEA HARENGUS",  # remove modeled species
                       "CLUPEIDAE"  # remove modeled species--assumed to be mostly atlantic herring
                       )) %>%
  #filter(!str_detect(pynam, "SHRIMP|CRAB")) %>%
  left_join(gencomlist) %>%
  filter(!gencom2 %in% c("ARTHROPODA", "ANNELIDA",
                         "CNIDARIA", "UROCHORDATA",
                         "ECHINODERMATA", "WORMS",
                         "BRACHIOPODA", "COMB JELLIES",
                         "BRYOZOA", "SPONGES",
                         "MISCELLANEOUS", "OTHER"))

fh.nefsc.pisc.pisccomplete.blueprey.unmod <- fh.nefsc.pisc.pisccomplete %>%
  mutate(blueprey = case_when(pynam %in% blueprey.unmod$pynam ~ "blueprey",
                              TRUE ~ "othprey"))

preystn.pisccomplete.unmod <- fh.nefsc.pisc.pisccomplete.blueprey.unmod %>%
  group_by(year, season, station) %>%
  count(blueprey) %>%
  pivot_wider(names_from = blueprey, values_from = n) %>%
  filter(year>1984)

#dim(preystn)[1]

bluepreystn.pisccomplete.unmod <- preystn.pisccomplete.unmod %>% 
  arrange(desc(blueprey)) %>%
  filter(!is.na(blueprey))

#dim(bluepreystn)[1]

Assign station id, change the months for spring and fall to align with the assessment, and calculate mean bluefish prey (for unmodeled species) per station in NEFSC:

bluepyall_stn.unmod <- fh.nefsc.pisc.pisccomplete.blueprey.unmod %>%
  #create id linking cruise6_station
  #create season_ng spring and fall Spring=Jan-May, Fall=June-Dec
  mutate(id = paste0(cruise6, "_", station),
         year = as.numeric(year),
         month = as.numeric(month),
         season_ng = case_when(month <= 6 ~ "SPRING",
                               month >= 7 ~ "FALL",
                               TRUE ~ as.character(NA))
         ) %>%
  select(year, season_ng, id, stratum,
         pynam, pyamtw, pywgti, pyvoli, blueprey, 
         pdcomnam, pdid, pdlen, pdsvol, pdswgt, 
         beglat, beglon, declat, declon, 
         bottemp, surftemp, setdepth) %>%
  group_by(id) %>%
  #mean blueprey g per stomach per tow: sum all blueprey g/n stomachs in tow
  mutate(bluepywt = case_when(blueprey == "blueprey" ~ pyamtw,
                               TRUE ~ 0.0),
         bluepynam = case_when(blueprey == "blueprey" ~ pynam,
                               TRUE ~ NA_character_)) 

# save at prey disaggregated stage for paper
saveRDS(bluepyall_stn.unmod, here("fhdat/bluepyall_stn_unmod.rds"))

stndat.unmod <- bluepyall_stn.unmod %>%
  select(year, season_ng, id, 
         beglat, beglon, declat, declon, 
         bottemp, surftemp, setdepth) %>%
  distinct()

#pisc stomachs in tow count pdid for each pred and sum
piscstom.unmod <- bluepyall_stn.unmod %>%
  group_by(id, pdcomnam) %>%
  summarise(nstompd = n_distinct(pdid)) %>%
  group_by(id) %>%
  summarise(nstomtot = sum(nstompd))

#mean and var pred length per tow
pisclen.unmod <- bluepyall_stn.unmod %>%
  summarise(meanpisclen = mean(pdlen),
            varpisclen = var(pdlen))

bluepyagg_stn.unmod <- bluepyall_stn.unmod %>%
  summarise(sumbluepywt = sum(bluepywt),
            nbluepysp = n_distinct(bluepynam, na.rm = T),
            npreysp = n_distinct(pynam),
            npiscsp = n_distinct(pdcomnam)) %>%
  left_join(piscstom.unmod) %>%
  mutate(meanbluepywt = sumbluepywt/nstomtot) %>%
  left_join(pisclen.unmod) %>%
  left_join(stndat.unmod)

# save at same stage as before, writing over old file
saveRDS(bluepyagg_stn.unmod, here("fhdat/bluepyagg_stn_unmod.rds"))

# current dataset, fix declon, add vessel
nefsc_bluepyagg_stn.unmod <- readRDS(here("fhdat/bluepyagg_stn_unmod.rds")) %>%
  mutate(declon = -declon,
         vessel = case_when(year<2009 ~ "AL",
                            year>=2009 ~ "HB", 
                            TRUE ~ as.character(NA)))

Combine with NEAMAP

Jim Gartland checked the NEAMAP dataset and found that silver hake and Atlantic mackerel are already absent from the prey list. Clupeidae may contain some small portion of Atlantic herring, but also includes river herrings and shads which we would want to keep in. So for NEAMAP we will keep the input data as is including Clupeids, but in the offshore dataset we assume Clupeids are mostly Atlantic herring so we leave them out.

NEAMAP inputs were similarly updated based on this predator size category list. NEAMAP includes the following predators, adding two not captured by the NEFSC survey offshore and leaving out those from NEFSC not captured inshore:

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

Read in new inputs, align columns, merge datasets, correct single station with wrong data in original dataset, and save:

neamap_bluepreyagg_stn <- read_csv(here("fhdat/NEAMAP_Mean stomach weights_Bluefish PreyWQ2.csv")) %>%
  mutate(vessel = "NEAMAP") %>%
  rename(id = station,
         sumbluepywt = sumbluepreywt,
         nbluepysp = nbfpreyspp,
         #npreysp = ,
         npiscsp = npiscspp,
         nstomtot = nstomtot, 
         meanbluepywt = meanbluepreywt,
         meanpisclen = meanpisclen.simple, 
         #varpisclen = ,
         season_ng = season,
         declat  = lat,
         declon = lon,
         bottemp = bWT,
         #surftemp = , 
         setdepth = depthm) 

  
# combine  
bluepyagg_stn.unmod <-  nefsc_bluepyagg_stn.unmod %>%
  bind_rows(neamap_bluepreyagg_stn) 

# check for incorrect NEAMAP station
bluepyagg_stn.unmod %>% filter(id == "NM20070901011") # has this station
# if sumbluepywt is 106564.2, this is incorrect
# corrected by Jim Gartland in September 2022

# correct single NEAMAP station 
bluepyagg_stn.unmod$sumbluepywt[bluepyagg_stn.unmod$id == "NM20070901011"] <- 4.8404
bluepyagg_stn.unmod$meanbluepywt[bluepyagg_stn.unmod$id == "NM20070901011"] <- 0.186169231


saveRDS(bluepyagg_stn.unmod, here("fhdat/bluepyagg_stn_all_unmod.rds"))

Merge with SST data to fill gaps

Now we take this dataset and add back the OISST data for stations missing surface temperature to apply our catchability covariates in VAST.

We can join this new dataset with the OISST dataset made previoulsy for survey stations, because we have not changed any station locations, only what we classified as bluefish prey at each station.

#read in diet data and station-OISST

bluepyagg_stn_all.unmod <- readRDS(here("fhdat/bluepyagg_stn_all_unmod.rds"))
dietstn_OISST <- readRDS(here("data-raw/dietstn_OISST.rds"))

dietstn_OISST_merge <- dietstn_OISST %>%
  dplyr::rename(declon = declon.x,
         declat = declat.x,
         year = year.x,
         oisst = sst) %>%
  dplyr::select(id, oisst) %>%
  sf::st_drop_geometry()

bluepyagg_stn_all_OISST.unmod <- left_join(bluepyagg_stn_all.unmod, dietstn_OISST_merge)

saveRDS(bluepyagg_stn_all_OISST.unmod, here("fhdat/bluepyagg_stn_all_OISST_unmod.rds"))

References

Garrison, L., and Link, J. 2000. Dietary guild structure of the fish community in the Northeast United States continental shelf ecosystem. Marine Ecology Progress Series, 202: 231–240. http://www.int-res.com/abstracts/meps/v202/p231-240/ (Accessed 22 October 2018).