The mesopelagic zone remains one of the most under-explored regions of the ocean; information about these deep-sea communities typically drops off with depth, despite their global distribution (Webb et al., 2010;St. John et al., 2016). Recent discoveries in the mesopelagic (200 -1,000 m depth) have provided new knowledge of pelagic invertebrate and fish communities (Kaartvedt et al., 2019) and resources for protein and nutraceuticals (St. John et al., 2016). Discovery of new taxa, ranging from vertebrates to viruses, continues (Robison, 2009;St. John et al., 2016;Govindarajan et al., 2021).
Mesopelagic fishes are found in all major oceans and are key players in food web dynamics of pelagic ecosystems (Williams et al., 2001;Kelly et al., 2019;Iglesias et al., 2023). Questions remain about biodiversity and biomass of mesopelagic fishes (Irigoien et al., 2014;Kaartvedt et al., 2019;Proud et al., 2019;Pauly et al., 2021) and their role in pelagic food webs and carbon cycling (Goetsch et al., 2018;McMonagle et al., 2023;Iglesias et al., 2023). Some species feed on particulate organic matter generated near the surface that sinks into the depths. Some epipelagic predators feed on migrating mesopelagic fishes at night by waiting for diel vertical migration (DVM) to bring their prey to them, while others dive to mesopelagic depths during the day to feed on these mesopelagic prey (Choy et al., 2017;Robison et al., 2020;Arostegui et al., 2022;Braun et al., 2022). Mesopelagic invertebrates and fishes mediate the relationship between primary production in surface waters, primary consumers in the epipelagic community, and deep-sea productivity and biomass; they are essential for the health and survival of deep-sea ecosystems (Choy et al., 2012(Choy et al., , 2013;;Irigoien et al., 2014;Käse et al., 2021). A more complete understanding of trophic dynamics of mesopelagic communities is critical considering anticipated growth in commercial harvesting of fish and other resources in the deep sea (St. John et al., 2016;Kelly et al., 2019), as well as possible impacts of climate change, plastic debris, deep-sea mining, and oil spills (Robison, 2009;Gamfeldt et al., 2015;Lusher et al., 2016;Levin et al., 2020;Kourantidou and Jin, 2022;Morzaria-Luna et al., 2022).
The diversity and abundance of mesopelagic fishes is exceptional, with species-specific variability in vertical distribution, DVM behavior, and trophic relationships (Caiger et al., 2021). This study examines seven species of fish collected in the Slope Water of the Northwest Atlantic Ocean (Figure 1). Argyropelecus aculeatus, silver hatchetfish, belongs to one of the most diverse, abundant, and widespread groups of mesopelagic fishes. The species has been reported to feed at dusk and to migrate into the epipelagic zone at night to feed (Hopkins and Baird, 1985). They are active predators and are known to feed on zooplankton, including copepods, ostracods, euphausiids, and gastropods, as well as gelatinous taxa and fish (Eduardo et al., 2020a;Receveur et al., 2020). The species has a typical life span of up to two years, with an average size of 58 mm; specimens up to 82 mm have been observed (Eduardo et al., 2019). Chauliodus sloani, Sloan's viperfish, is a predator equipped with massive teeth to cage trapped prey (Bergman et al., 2023). This species occurs in oceans throughout the world and has been reported to exhibit DVM. Reported depth ranges are variable, perhaps due to differences in behaviors among life history stages, including daytime depths below 400 m (Eduardo et al., 2019(Eduardo et al., , 2020b)), with night-time migration to within 100 m of the surface (Sutton and Hopkins, 1996) to feed on smaller prey (Butler et al., 2001). The species has been found to prey mainly on fishes up to half of their body size, with larger prey captured by the unhinged jaw and digested by an elastic stomach (Battaglia et al., 2018). Hygophum hygomii, Bermuda lanternfish, also shows DVM behavior and is known to feed on zooplankton, including primarily copepods, amphipods, and euphausiids (Pusch et al., 2004). Typical size ranges are 29-63 mm, with maximum size of 70 mm (Pusch et al., 2004). Malacosteus niger, stoplight loosejaw, is highly adapted to the deep-sea environment, with a large gape and massive fangs. The species has been observed to be up to 270 mm in length and is found worldwide. The species has not been observed to exhibit DVM. Known prey include zooplankton, especially copepods, and larger micronekton, including decapods and other fishes (Sutton, 2005). Nemichthys scolopaceus, slender snipe eel, occurs throughout temperate to tropical regions of all major oceans. The species is distributed between 200-1000 m and exhibits vertical migration patterns, although the behavior may be caused by following prey (Smith and Tighe, 2002;Mundy, 2005). Prey are mostly large crustaceans, with occasional smaller zooplankton, pelagic mollusks, and fish. The species has been observed reaching lengths up to 1150 mm (Fraser-Brunner, 1936;Feagans-Bartow and Sutton, 2014). Scopelogadus beanii, Bean's bigscale, does not exhibit DVM and feeds at depth; the diet is comprised of micronekton, including amphipods and gelatinous zooplankton, and fish (Gartner and Musick, 1989). Sigmops elongatus, elongated bristlemouth, is known to undergo DVM and prey on fishes and larger invertebrates; the species has been recorded up to 280 mm in length (Woodstock et al., 2020).
Integrative molecular (metabarcoding) and morphological (microscopic) analysis of fish diets has been identified as a best practice in several studies (Amundsen and Sańchez-Hernańdez, 2019;Clarke et al., 2020). Recommended approaches for quantitative morphological analysis of prey include comparisons of counts, frequencies of occurrence, or inferred biomass for selected prey groups (Matthews et al., 2021), using methods developed for integrative analysis of pelagic biodiversity (Ershova et al., 2021). Metabarcoding is defined as the identification of multiple taxonomic groups or species from a bulk sample containing genetic material from diverse sources (Taberlet et al., 2012). The emergence of metabarcoding was driven by the development and rapid advances in high-throughput DNA sequencing (HTS), which can yield billions of sequences and rapidly analyze unsorted samples from any environment or ecosystem (Taberlet et al., 2012). The marine metabarcoding revolution has been built upon several global initiatives, including characterization of biodiversity of pelagic ecosystems using ribosomal RNA (rRNA) genes (Quast et al., 2013) and species identification of marine metazoa using the mitochondrial cytochrome oxidase I (COI) barcode region (Bucklin et al., 2011(Bucklin et al., , 2021b)). A necessary foundation and key factor for applications of metabarcoding is the availability of reference databases to enable taxonomic identification of specimens by matching metabarcode sequences to reference sequences from specimens identified by morphological taxonomic experts (Coissac et al., 2012;Alberdi et al., 2017). Concerted efforts are continuing toward taxonomically-complete and geographically-comprehensive reference sequence databases for marine organisms, including zooplankton and fish (Sato et al., 2018;Blanco-Bercial, 2020;Bucklin et al., 2021b;Govindarajan et al., 2023).
Metabarcoding has been used to examine the diets of fishes from diverse ocean habitats, including Acanthopagrus latus, yellowfin seabream, in coastal waters of China (Su et al., 2018;Pan et al., 2021), Boreogadus saida, polar cod, in the Barents Sea (Maes et al., 2022), Syngnathus watermeyeri, and S. temminckii, estuarine pipefishes in South Africa (Serite et al., 2023), Lethenteron camtschaticum, Arctic lamprey, in the Bering Sea (Shink et al., 2019), Sardina pilchardus, European sardine, and Sprattus, European sprat, in the Bay of Biscay (Albaina et al., 2016), and numerous species of coral reef fishes (Casey et al., 2019). The diets of mesopelagic fishes have been analyzed using metabarcoding in Antarctic waters by Clarke et al. (2020). Recent efforts have explored broader use of metabarcoding of fish stomach contents for applications in fisheries and ecosystem management (Canals et al., 2024).
Detection and identification through morphological (microscopic) analysis of prey species in the gut contents of fish has been widely used to examine trophic relationships of marine and freshwater species (Baker et al., 2014;Young et al., 2015;Suca et al., 2018;da Silveira et al., 2020). Comprehensive reviews have considered the impacts and implications of diverse analytical approaches (Manko, 2016;Amundsen and Sańchez-Hernańdez, 2019). Recent review papers have recommended analytical or statistical approaches (Amundsen and Sańchez-Hernańdez, 2019) and standardization of methodologies (Buckland et al., 2017). Questions remain about best analytical approaches, including counts of numerical abundance, frequency of occurrence, and volume or biomass, as well as the most appropriate statistical tests of the resulting data, which include various indices of relative importance and selectivity (Manko, 2016).
This study characterizes and compares the diets of seven cooccurring mesopelagic fish species in the NW Atlantic Slope Water using integrative analysis. Our analyses explore the power of metabarcoding for detection and identification of invertebrate prey species. Morphological (microscopic) counts of prey in gut contents of fish collected in the same net samples allowed further consideration of the benefits and challenges of integrative analysis for understanding the trophic dynamics of the mesopelagic zone.
Fish samples were obtained from the Slope Water of the Northwest Atlantic Ocean in August, 2018 and July-August, 2019 during cruises of the NOAA research vessel (R/V) Henry B. Bigelow (Figure 2A). Hydrographic data were collected during both cruises, with multiple conductivity, temperature, depth (CTD) profiles at each station (Figure 2B). Sampling for invertebrates and fish was done using a modified Marinovich midwater trawl (De Robertis et al., 2017;Jech and Lavery, 2018) and a 1-m 2 MOCNESS (Wiebe et al., 1985), with deployments to differing depths and times of day or night (Table 1). Specimens planned for analysis of gut contents Mesopelagic fish species analyzed using metabarcoding of DNA extracted from dissected gut contents. Specimens were collected in 2018 and 2019; see Table 1 for collection information and Table 4 for sample sizes. Photos by Paul Caiger (University of Auckland).
were identified to species, measured (standard length, mm), immediately flash-frozen in liquid nitrogen to stop digestive processes, and stored at -80°C. Up to 10 fish per species per trawl were flash frozen.
Post-cruise processing of the frozen fish from each midwater trawl was carried out in a research laboratory at the Woods Hole Oceanographic Institution (WHOI, Woods Hole, MA). Fish were photographed, measured for standard length, fork length, and total length (in mm), weighed, and then dissected. Muscle tissue was excised from representative specimens to identify the species based on DNA sequencing of the COI barcode gene region for comparison to reference databases (Ward et al., 2009;Steinke and Hanner, 2011;Govindarajan et al., 2023;Quigley et al., 2023). Specimens were partially thawed to allow dissection of guts; gut contents were removed and weighed and transferred to individual vials. Gut contents were either preserved in 95% ethanol for morphological identification or frozen at -80°C for metabarcoding. During both HB-2018 and HB-2019 cruises, every attempt was made to allocate fish from each net tow equally for DNA metabarcoding and morphological microscopy based upon both numbers and sizes of specimens (Supplementary Tables 2, 3).
Gut contents of each fish were thawed on ice in elution buffer made of Tris-HCl (10mM), EDTA (100 mM, pH 8), NaCl (200 mM), sodium dodecyl sulfate (SDS, 1%), and Milli-Q water. After homogenization for 15 sec, samples were incubated in a water bath at 55°C for 6-7 hours. DNA was extracted from digested samples using phenol:chloroform:isoamyl alcohol (25:24:1) and cleaned using a DNEasy PowerClean Pro Kit (Qiagen) according to manufacturer instructions.
Purified DNA from gut contents was used to amplify two gene regions: V9 hypervariable region of 18S rRNA (Amaral-Zetter et al., 2009) and a portion of COI (Leray et al., 2013). All forward and reverse primers were altered for multiplexed sequencing by adding 5' adapters (Illumina, Inc., San Diego, CA). The V9 18S rRNA reaction used 4 µL of DNA template, with KAPA HiFi reagents (KAPA Biosystems, Massachusetts, USA): 5 µL buffer containing MgCl 2 , 1µL dNTPs, 0.5 µL HiFi Taq Polymerase, and 1 µL of each primer (10 µM). The primer pair 1380F and 1510R (Amaral-Zetter et al., 2009) was used to amplify V9 18S rRNA, with the following PCR protocol: one denaturation cycle at 98°C for 30 sec; 10 cycles of 98°C for 20 sec, 56°C for 30 sec, and 72°C for 15 sec; and 16 cycles of 9°C for 10 sec, 66°C for 30 sec, and 72°C for 15 sec; and 1 extension cycle of 72°C for 7 min. PCR amplification of a 313 base-pair (bp) region of COI used the primer pair mlCOIintF and jgHCO2198 (Geller et al., 2013;Leray et al., 2013). The PCR reaction used 20 ng of DNA, with Platinum Taq reagents, 4 µL buffer, 2.4 µL MgCl 2 , 0.8 µL dNTPs, 0.2 µL HiFi Taq Polymerase, and 0.8 µL of each primer (10 µM), with the following protocol: one denaturation cycle at 94°C for 60 sec; 38 cycles of 94°C for 30 sec, 46°C for 30 sec, 72°C for 90 sec; and 1 extension cycle of 72°C for 5 min. Both V9 18S rRNA and COI amplicons were checked for successful amplification by running in a 2% agarose gel with a 50 base-pair (bp) marker.
Library preparation entailed adding index primers in a second PCR amplification of the purified amplicons using a master mix composed of (per sample): 5.0 µL purified PCR product; 5 µL Nextera XT Index 1 Primer; 5 µL Nextera XT Index 2 Primer; 25 µL 2x KAPA HiFi HotStart ReadyMix; 10 µL PCR-grade water; for a total volume of 50 µL. The PCR protocol was: 95°C for 3 min; 8 cycles of: 95°C for 30 sec, 55°C for 30 sec, 72°C for 30 sec; and 1 cycle of 72°C for 5 min. The indexed PCR product was purified using AMPure XP beads, with a final elution volume of 25 µL. Successful library attachment was verified using an Agilent 2200 TapeStation automated electrophoresis system. Libraries were quantified using a Qubit 3.0 fluorometer, normalized according to amplicon size, pooled, and denatured with 0.2 N NaOH. Samples were spiked with a minimum of 5% PhiX (Illumina, Inc.). Bidirectional sequencing was carried out at the University of Connecticut Center for Genomic Innovation (CGI; https:// cgi.uconn.edu/) on an Illumina MiSeq sequencer using the MiSeq Reagent Nano Kit Ver. 2 (500 cycles; 1 million clusters) spiked with a minimum of 5% PhiX.
Demultiplexed reads for V9 were processed using a custom script (URL available upon publication) for the Mothur pipeline (Ver. 1.44.3;Schloss et al., 2009) and run on the Xanadu computing cluster of the UConn Computational Biology Core (CBC; https:// bioinformatics.uconn.edu/). Contiguous sequences (contigs) were assembled from forward and reverse Illumina MiSeq reads and trimmed to the overlapping section. Sequences containing ambiguous bases, quality Phred scores < 30, and with lengths shorter than 120 bp were removed from analysis. Unique sequences were aligned against the reference database, SILVA Release 132 (Quast et al., 2013; https://www.arb-silva.de/ documentation/release-132/). Sequences were trimmed to a uniform length by removing the beginning and terminal ends of sequences that extended beyond the V9 gene region. Any sequences that did not span the entire V9 region were removed, decreasing the likelihood of artefactual operational taxonomic units (OTUs) being created during clustering. Concerns that PCR error may contribute to errors in biodiversity assessment (Kelly et al., 2019) were addressed by using the UNOISE 37 method (Edgar, 2016) within Mothur (Ver. 1.44.3) to de-noise aligned sequences before clustering of OTUs, which was done using a limit of 2 bp difference between sequences. Sequences were screened for chimeras using the VSEARCH command (Rognes et al., 2016); sequences with chimeras were removed from analysis.
OTUs were assigned taxonomic identifications using a tailored 18S rRNA database by Blanco-Bercial (2020) that adds sequences for eukaryotic marine organisms acquired from the NCBI GenBank sequence repository to the SILVA 132 database (Quast et al., 2013). Taxonomic assignments were determined using a naïve Bayesian classifier algorithm (Wang et al., 2007), which uses the highest probability that a given sequence contains kmers (i.e., DNA segments of length k in nucleotides) specific to a sequence of a known taxonomic identity; default kmer size (ksize) = 8 was used. Taxonomic assignments were based on bootstrap values ≥ 80% after 100 iterations.
Bi-directional COI reads were processed using a custom script for Mothur (Ver. 1.44.3;Schloss et al., 2009) and run on the Xanadu computing cluster. The bioinformatics pipeline used to quality control and filter COI sequences closely followed the 18S rRNA Mothur script. Contigs greater than 150 bp in length were aligned to the Global MetaZooGene Database (MZGdb; https://metazoogene. org/MZGdb; accessed March 3, 2022), a custom COI database including publicly available sequences downloaded from GenBank and BOLD (Bucklin et al., 2021b). Sequences were further denoised of PCR and sequencing errors using the UNOISE3 method in Mothur (Edgar, 2016). Chimeras were identified using the VSEARCH command (Rognes et al., 2016) and removed from analysis.
Taxonomic identifications were based on sequences and OTUs for COI determined using a naïve Bayesian classifier algorithm in Mothur (Ver. 1.44.3) against the geographic-specific North Atlantic MZGdb (Bucklin et al., 2021b). Taxonomic assignments for specieslevel identifications used bootstrap values ≥ 97% after 100 iterations (Schroeder et al., 2020). Before performing zooplankton community analyses, sequences with abundances < 2 (i.e., global singletons) across the entire dataset were removed.
Specimens of seven mesopelagic fish species were analyzed by metabarcoding of DNA purified from gut contents. Total samples sizes were: 36 specimens for V9 and 28 specimens for COI from 2018 samples; 30 specimens for V9 and 33 specimens for COI in 2019 (Supplementary Table 1). Multivariate statistical analysis to characterize diet diversity within and among the mesopelagic fish species was carried out for sequence numbers and OTUs of both COI and V9. Results described here focus on sequence numbers. Sequence numbers are reported as average percent occurrence (Log10 + 1 transformed) for fish samples of each species collected in each of 2018 and 2019. Results are reported and analyzed as average proportions of sequence numbers of each gene region (Log10 + 1 transformed) for prey groups (major taxonomic categories) selected based on metabarcoding results as reported in the taxonomy summary file (Wang et al., 2007) generated by Mothur (Ver. 1.44.3;Edgar, 2016) for each analysis. Results were analyzed for eight taxonomic groups (Copepoda, Eucarida, Amphipoda, Ostracoda, Chaetognatha, Hydrozoa, Cephalopoda, and Gastropoda) for which sequences of both gene regions detected prey (i.e., non-zero results in at least one fish species). Analysis was done for an additional six prey groups (Polychaeta, Salpida, Nematoda, Platyhelminthes, Scyphozoa, and Ctenophora), which were detected only by V9. COI sequence numbers showed a higher frequency of zero values. Fish specimens showing zero sequences for either V9 or COI for all (or in several cases, almost all) prey groups were removed from analysis, reducing sample sizes in some cases.
V9 and COI sequences classified as fish (Teleostei) in gut contents were reported for each sample (Supplementary Tables 2A, B), but not included in the analyses, even for known piscivores, since neither gene region can discriminate sequences resulting from DNA of predator gut tissue versus ingested prey. COI sequence numbers for identified species of fish from gut content DNA were reported for each predator species (Table 2).
Ecological network analysis based on V9 was used to evaluate predatorprey linkages between the seven fish species and 14 prey groups (Table 3) using BitMAT (Flores et al., 2016). The percentages of V9 sequences for the fish predator species were computed for fish from both 2018 and 2019; prey with ≥ 1% were included in the analysis. Sequences were summed for each of the 14 prey groups, divided by the number of fish analyzed for each predator species, and multiplied by 100 to give percentages.
TABLE 2 Prey species found in fish guts based on COI sequence numbers.
Table 2A. Species of Crustacea detected by COI metabarcoding of fish gut contents. Fish Species Prey Species by Group Aa Cs Hh Mn Ns Sb Se Copepoda (N=22) Calanus finmarchicus --0.48 -0.33 11.29 -Calanus hyperboreus 0.21 -0.19 -0.17 --Calocalanus pavo --2.10 ----Candacia curta --1.29 ----Centropages hamatus 0.07 ------Clausocalanus lividus --0.14 ----Clausocalanus pergens --0.14 ----Euchirella splendens --1.29 ----Mesocalanus tenuicornis ----0.17 --Nannocalanus minor 0.07 -103.38 ---2.10 Oithona similis --0.10 ----Pareucalanus sewelli 0.14 -1.05 ----Pleuromamma abdominalis --1.86 ---1.30 Pleuromamma borealis 0.14 -3.52 -2.00 -0.80 Pseudocalanus moultoni --0.10 ----Scolecithrix danae 0.07 -192.10 -0.33 0.14 6.90 Spinocalanus antarcticus -----0.14 -Subeucalanus crassus -----10.71 -Temora longicornis --0.05 --0.14 -Temora stylifera --98.00 -0.33 --Temora turbinata 0.07 -0.10 ----Undinula vulgaris --1.57 ----Eucarida (N=18) Acanthephyra purpurea ----762.17 --Achelous ordwayi 18.57 -0.67 ----Belzebub faxoni --0.29 ----Calappa flammea --2.29 ----Callinectes sapidus --0.62 ----Euphausia americana 0.14 -63.81 ----Euphausia krohnii 92.79 35.71 13.00 ---1565.60 Euphausia tenera 0.07 -0.14 ----Eusergestes arcticus -0.18 -----Meganyctiphanes norvegica 274.29 74.29 0.05 ---617.80 Metapenaeopsis gerardoi -----0.14 -Nematoscelis megalops 1398.07 58.82 768.57 -0.17 0.14 1.00 Neosergestes edwardsii -10.94 -----(Continued) Table 2B. Species of non-crustacean groups detected by COI metabarcoding of fish gut contents. Fish Species Prey Species by Group Aa Cs Hh Mn Ns Sb Se Hydrozoa (N=4) -Obelia dichotoma 0.07 ------Lensia achilles -----6.86 -Nanomia cara -----0.14 -Rhopalonema velatum ---0.33 -2.86 -Cephalopoda (N=5) Haliphron atlanticus ---0.67 ---Histioteuthis reversa ---0.67 ---Ommastrephes bartramii --0.14 ----Pterygioteuthis gemmata --0.05 ----Sthenoteuthis pteropus 4.50 -0.05 --0.86 1680.20 Gastropoda (N=7) Atlanta selvagensis --0.05 ----Clio pyramidata 4.14 ---0.17 --Creseis conica --0.14 ----Diacavolinia longirostris --0.10 ----Diacria trispinosa 257.43 0.06 -----Firoloida desmarestia --2.57 ----Heliconoides inflatus 0.21 -0.33 --1.29 0.10 (Continued)
Patterns of variation in sequence numbers of COI and V9 for the selected prey groups were statistically evaluated among species and between years of collection in MatLab (Ver. 2020B). Two-way Analysis of Variance (ANOVA) was done to test for significant differences in sequence numbers (Log10 + 1 transformed) of both V9 and COI for selected prey groups between the predator fish species and the two years of collection (2018 and 2019). One distance measure used was Bray-Curtis dissimilarity coefficient (Bray and Curtis, 1957;McCune et al., 2002), with results displayed by cluster diagrams. Differentiation among the seven species was evaluated by Non-Metric Multidimensional Scaling (NMDS) using the FATHOM Toolbox for MatLab (Jones, 2017; https://www.usf.edu/marine-science/research/matlab-resources/ index.aspx/). The Shannon Diversity Index (H; Pielou, 1977) was calculated using numbers of sequences for the prey groups detected for each gene.
Identified prey species with a minimum threshold of 1,000 COI sequences in total across all predators were reported for 7 of the 8 prey groups analyzed for COI (Table 2), excluding Chaetognatha, for which no species were identified, and fish (Teleostei). Normalization was done separately for each s p e c i e s b a s e d o n o v e r a l l m e a n s ; s t a t i s t i c a l a n a l y s i s included NMDS.
Fish used for morphological (microscopic) analysis of gut contents were selected from the same collections and net tows during the 2018 and 2019 cruises of the R/V Henry B. Bigelow as the samples used for metabarcoding of gut contents (as described above). Gut contents for microscopy were thawed, placed in 95% ethanol for preservation, and identified to the lowest taxonomic unit practical, focusing on groups that could be accurately and reliably distinguished after partial digestion, using a Leica M60 dissecting microscope (Suca et al., 2018).
Morphological counts of prey in the gut contents of specimens of each fish species were classified into six of the taxonomic groups used for analysis of COI and V9 metabarcoding results: Copepoda, Eucarida, Amphipoda, Ostracoda, Chaetognatha, and Gastropoda. An additional group included prey fishes (Teleostei). A general category of soft-bodied organisms was used to record prey that could not be identified morphologically, including Chaetognatha, Hydrozoa, and Cephalopoda. Numbers of prey in each category were recorded for each individual gut content sample analyzed and averaged across all specimens examined for each predator species. Statistical analysis included calculation of the Bray-Curtis dissimilarity coefficient (Bray and Curtis, 1957;McCune et al., 2002); results were displayed by cluster diagrams. NMDS analysis used the FATHOM Toolbox for MatLab (Jones, 2017; https:// www.usf.edu/marine-science/research/matlab-resources/ index.aspx/). The Shannon Diversity Index (H; Pielou, 1977) was calculated using average counts for prey groups for each fish species.
The relationships between size of fish (standard length, mm) and Shannon Indices of prey diversity, based on COI and V9 metabarcoding data and morphological counts, were examined to determine whether fish size was a significant predictor of the diversity of prey ingested. Regression analyses were carried out in MatLab (Jones, 2017). Two fish species (A. aculeatus and H. hygomii) had enough data for regression analysis of size versus V9 and COI sequences and counts; one species was analyzed only for V9 (C. sloani); one species (S. elongatus) was analyzed for V9 and prey counts.
To examine relationships between prey diversity and time and depth of predator collection, CTD data were used to characterize temperature and salinity of the water column when and where the fish were collected in both sampling years. Temperature and salinity data from multiple CTD casts were pooled; maximum, minimum, and average values were determined for selected depth ranges and graphed on the same plot to allow comparisons between CTD casts and years. Variation of V9 and COI sequence numbers and morphological prey counts based on both average values for the seven fish species and individual fish specimens were examined in relation to depth of the net tow and hydrographic structure of the water column at the time of collection. Impacts of collection time on prey composition focused on comparison of samples from day and night tows during 2019 (Table 1).
Hydrographic conditions at the sampling sites in the NW Atlantic Slope Water were described based on multiple CTD profiles in both 2018 and 2019 (Figure 2B). Temperatures and salinities at the sample collection locations in HB-1805 (Station #2) and HB-1907 (Stations #2 and #3) were similar below 600 m. At mid-depths (50-600 m), both temperatures and salinities were higher and more variable at Station 3 in 2019 than at Station 2 in either year. Near-surface temperatures (above 50 m) were much warmer at all three stations in both 2018 and 2019, but relatively cooler in 2018 at Station 2. Above 50 m, salinities showed wide variation at Station 2 in both years, but were lower in 2018, with higher and less variable values in 2019.
V9 metabarcoding detected the presence of 14 prey groups in the gut contents of the seven fish predators (Table 3; Supplementary Table 2A). V9 sequences for fish prey were excluded from analysis The relative importance of the prey groups for each predator species is shown in a bipartite graph (Dormann et al., 2009) indicating quantitative linkages between predators and prey (Figure 3). The network diagram clearly shows the importance of three prey groups, Eucarida, Amphipoda, and Copepoda, which have strong links, shown as thick bars, to five predators: A. aculeatus, H. hygomii, M. niger, and S. elongatus. Another prey group, Salpida, has links of varying width to all predators except H. hygomii (Figure 3).
Statistical analysis of 868 individual observations by 2-way ANOVA revealed significant variation among the predator
Network diagram of frequencies of 14 prey groups based on V9 18S rRNA sequences for 7 fish predator species. The width of each link is proportional to the frequency of prey group occurrences in each fish species, computed by adding sequences for individual fish samples, dividing by the number of fish analyzed from 2018 and 2019, and multiplying by 100. Links to prey with frequencies < 1% are not shown (i.e., Ostracoda, Chaetognatha, Scyphozoa). See Table 3 for explanation of abbreviations for predator species.
species (p < 0.001) and between years (p < 0.043), with significant interaction (p < 0.002). Bray-Curtis dissimilarity cluster plots based on V9 for prey groups (Figure 4A) showed differences among the fish species, with clustering of H. hygomii from 2018 and 2019, along with M. niger, and separate-but-overlapping clusters of A. aculeatus and C. sloani; S. elongatus clustered by year, while M. niger, N. scolopaceus, and S. beanii showed lack of separation. Nonmetric multidimensional scaling (NMDS) based on V9 detection of prey groups for each year also revealed variation among the seven fish species, with distinct groupings of H. hygomii and S. beanii in 2018, with small and distinct clusters of A. aculeatus and H. hygomii, and a broad group for C. sloani in 2019 (Figure 5). The Shannon Diversity Index based on V9 showed variation in prey groups between predators for both 2018 and 2019, albeit with significant variation among fish in both 2018 and 2019 (Figure 6A).
COI metabarcoding yielded sequences that were classified to eight prey groups (Table 4; Figure 7); other groups detected by V9 showed zero COI sequences for all fish analyzed. Fish were not included in analysis; COI sequence numbers for fish as a prey group were very large for all predators (Supplementary Table 2B). Gut contents of A. aculeatus showed COI sequences for 7 of the 8 prey groups detected in at least one sample over the two years, albeit mostly in low sequence numbers. The prey group showing highest COI sequences was eucarids, which were abundant in fish collected during night-time tows in 2018 and 2019. Of 17 C. sloani from both years, COI detected prey in only one specimen, with eucarid sequences from one sample collected during a 2018 night-time tow. Hygophum hygomii showed high COI sequence numbers for Bray-Curtis similarity cluster plots for fish gut contents of samples collected during 2018 and 2019 based on (A) V9 sequence numbers (Log10 + 1) for 14 prey groups, (B) COI sequence numbers (Log10 + 1) for 8 prey groups, (C) microscopic counts for prey groups, soft-bodied organisms and fish. Sample labels indicate species and year of collection, either 2018 (18) or 2019 (19). Bucklin et al. 10.3389/fmars.2024.1411996 Frontiers in Marine Science frontiersin.org eucarids and copepods in most samples from both 2018 and 2019; amphipods were also abundant in three fish from 2018 night-time tows and one fish from a 2019 night-time tow. Small COI sequence numbers were found in three M. niger collected in a 2018 night-time tow, including eucarids, amphipods, hydrozoans, and cephalopods; no fish were analyzed for 2019. COI detected low levels of copepods and eucarids in gut contents of N. scolopaceus for both years, with high numbers for one fish collected in a 2019 day-time tow. Most of seven S. beanii had very low or zero COI sequence numbers for most prey groups, with the exception of one sample collected during a 2019 night-time tow, which included copepods, ostracods, and hydrozoans. Eucarids were detected in gut contents of four of 10 S. elongatus; copepods were found in one fish from a 2018 night-time tow and cephalopods in one fish from a 2019 night-time tow.
The two-way ANOVA based on 237 individual observations for COI sequences for four predators showed no significant variation among species (p < 0.710) or years of collection (p < 0.632). Three predator species were not included in this analysis: M. niger was not caught in 2019; only one fish for each of C. sloani and N. scolopaceus had non-zero COI sequence numbers for the eight prey groups.
Patterns of diet variation among the predators revealed by Bray-Curtis clustering based on COI sequence numbers showed differences among species and years, but with no clear clustering patterns (Figure 4B). Most notable was the similarity of four H. hygomii collected in 2019 and clustering of A. aculeatus from both 2018 and 2019. NMDS analysis showed no clear patterns within and among predator diets for 2018 and 2019 based on COI, with the exception that H. hygomii clustered together in 2019 (Figure 8). Shannon Diversity Index for fish diets ranged widely for most species, with the exceptions of M. niger, based on only two fish, and C. sloani and S. elongatus, which both showed low diversity (Figure 6B).
COI metabarcoding of predator gut contents identified a total of 77 prey species, including 49 species of crustaceans (Table 2A). Other prey groups for which species were identified by COI included hydrozoans, cephalopods, gastropods, and fish (11 species, not including predators; Table 2B). Three euphausiid species, Euphausia krohnii, Meganyctiphanes norvegica, and Nematoscelis megalops, were important prey for A. aculeatus, C. sloani, H. hygomii, and S. elongatus. The predator with the highest number of prey species detected was H. hygomii, with a total of 46 species, including 18 species of copepods, 12 species of eucarids; and seven species of amphipods. COI identified several fish species (excluding predators) as prey of M. niger, including Benthosema glaciale and Notoscopelus resplendens. No M. niger were analyzed for 2019. For gut contents of S. beanii, COI detected the copepod species, Calanus finmarchicus and Subeucalanus crassus in the highest frequencies; an ostracod, Conchoecissa ametra, occurred in one sample from 2019.
The highest total average sequence numbers for prey species across all predators were the euphausiids, Euphausia krohnii, Meganyctiphanes norvegica, and Nematoscelis megalops; and the squid, Sthenoteuthis pteropus (Table 5A). For all predator fish samples analyzed from 2018 and 2019 collections, the most usual finding was zero sequences for species of hydrozoans, cephalopods, and gastropods. The highest numbers of COI sequences for prey species in any one predator were the fish, Paralepis coregonoides and Scopelosaurus lepidus, found in A. aculeatus. The COI metabarcoding results included exceptional numbers of sequences for the predators: for six of the seven species, by far the largest number of COI sequences for any species (>10,000) matched the predator (Table 2B). COI cannot distinguish consumed prey from the fish gut tissue (i.e., predator swamping; Albaina et al., 2016;Pereira et al., 2017;Clarke et al., 2020). No correlation was found between COI sequence numbers for fish species and numbers in each sample, suggesting that neither "codend feeding" nor damage during collection were occurring, although prey of prey and eDNA sources cannot be excluded.
Morphological counts of prey were summarized by taxonomic group and averaged across all specimens examined for gut contents Non-Metric Multi-Dimensional Scaling (NMDS) analysis of V9 sequence numbers (Log10 + 1) for 14 prey groups in fish gut contents for samples collected in 2018 and 2019. See Table 4 for list of prey groups included in the analyses.
of each predator (Table 5; Supplementary Table 3). Counts of A. aculeatus found euphausiids, amphipods, gastropods, and hydrozoans in many of the specimens examined, with highest frequencies of amphipods and soft-bodied prey in 2018 and gastropods in 2019. Counts of C. sloani detected a copepod and soft-bodied prey in two specimens collected in 2018. Most samples of C. sloani examined showed empty guts based on morphological counts. Counts for H. hygomii showed highest averages for copepods and amphipods in 2018 (Table 5A); in 2019, copepods were detected in all H. hygomii examined (Supplementary Table 3). Crustaceans predominated in M. niger, with eucarids detected in both years and copepods in 2019 (Table 5). Microscopic analysis of N. scolopaceus found almost no prey, except for one amphipod in one sample. Crustacean groups and soft-bodied prey were detected in S. beanii each year, with the exception of eucarids in 2019. Morphological counts for S. elongatus reported the same crustacean groups (excluding ostracods) and soft-bodied prey in both years.
Based on proportional numbers of prey groups discriminated by microscopy, there were some clear differences among the predators (Figure 9). In both 2018 and 2019, soft-bodied prey were found frequently in A. aculeatus, S. beanii, and S. elongatus, rarely in H. hygomii, and were absent from M. niger and N. scolopaceus (Table 5; Figure 9). Statistical analysis of counts by species and years was not possible by ANOVA due to many zero values. The Bray-Curtis coefficient (Bray and Curtis, 1957;McCune et al., 2002) and resulting cluster diagram showed intermixing of diet diversity for all 7 species (Figure 4C). The NMDS analysis showed similarity of prey composition for A. aculeatus and H. hygomii in both 2018 and 2019, albeit with some outliers (Figure 10). The Shannon Diversity Index (Pielou, 1977) showed broad ranges for species with the largest sample sizes, A. aculeatus and H. hygomii (Figure 6C).
Statistical tests of the relationship between sizes of fish sampled (measured as standard length, SL) and the Shannon Index of prey diversity based on V9, COI, and morphological counts revealed no significant relationships based on multi-dimensional regression analysis. None of the regressions were significant, indicating no relationship between the diversity of prey and fish standard length for any predator collected over the two years. The similarity of the hydrographic structure at the collection sites in 2018 and 2019 (Figure 2) allowed direct comparisons between years. Statistical analyses carried out to evaluate variation in prey group composition based on V9 and COI in relation to collection parameters, including time (day versus night) and depth of tow, did not reveal any significant relationships. Patterns of variation of prey species diversity based on COI also showed no consistent relationships with time of collection or depth of the tow for any predator.
Integrative analysis of mesopelagic fish diets Mesopelagic fishes are key players in deep sea food webs and have significant impacts on ecosystem function. They are known to be important prey for higher trophic levels (Iglesias et al., 2023) and transporters of organic material in the open ocean (Saba et al., 2021), although knowledge is incomplete of their trophic relationships and the composition of their diets (Goetsch et al., 2018;Clarke et al., 2020;Robison et al., 2020). Molecular approaches are widely recognized to have yielded new understanding of trophic interactions in ocean ecosystems (Symondson and Harwood, 2014), with promise for continued successful development (Roslin et al., 2019). An important contribution of metabarcoding for analysis of the diets of mesopelagic fishes is the confirmation of the importance of gelatinous zooplankton, which are difficult to detectand certainly identifyin morphological analysis of fish gut contents (Choy et al., 2017). V9 can provide accurate and semi-quantitative (i.e., proportional frequencies) detection of prey across a broad phylogenetic spectrum, including microbes, protists, and metazoans (Bucklin et al., 2019;Blanco-Bercial, 2020;Govindarajan et al., 2021). Hydrozoans and salps were shown to be significant food sources for all predators, while ctenophores were abundant in several species (A. aculeatus, C. sloani, S. beanii). COI also revealed species of somebut not allgelatinous groups, including three species of hydrozoans, Lensia achilles, Nanomia cara, and Rhopalonema velatum, in S. beanii; Sthenoteuthis pteropus (Cephalopoda) was identified as a primary prey of S. elongatus and also found in A. aculeatus. These findings confirmed the broad prey diversity of some mesopelagic fish species and provided detailed information supporting the importance of gelatinous zooplankton in mesopelagic food webs (Purcell and Sturdevant, 2001;Choy et al., 2017).
Metabarcoding analyses that identify prey species can also support the important goal of determining sources of primary productivity supporting mesopelagic biomass. For instance, the primary prey species found in gut contents of H. hygomii included the surface-dwelling, non-migrating copepod, Nannocalanus minor (Ambler and Miller, 1987) and the middepth, non-migrating euphausiid, Nematoscelis megalops (Endo and Wiebe, 2007). Although H. hygomii exhibits DVM and the net tows varied in both time and depth, there were no apparent trends in prey for day versus night. Our findings suggest that this species feeds predominantly in surface waters at night with little feeding at their daytime depth. Another predator, A. aculeatus, showed large numbers of COI sequences for two euphausiids: nonmigrating N. megalops and migrating Meganyctiphanes norvegica (Tarling et al., 1999). Large numbers of COI sequences were found in N. scolopaceus for both migrating and non-migrating prey: the decapod shrimp, Acanephyra purpurea, which shows DVM (Roe, 1984) and a non-migrating deep-dwelling fish, Paralepis coregonoides (Harry, 1953). The gut contents of the nonmigrating predator S. beani contained representatives of all prey groups analyzed, including both migrating and non-migrating species with extensive (>1,000 m) vertical ranges: the hydrozoan Rhopalonema velatum (Lucǐćet al., 2009); ostracod Conchoecissa ametra (Chavtur and Mazdygan, 2015); copepod Calanus finmarchicus (Fiksen and Carlotti, 1998); and fish Scopelosaurus lepidus (Gordon and Duncan, 1985). Taken together, our data revealed a food web that is apparently well-connected across depth zones in the open ocean (Sutton, 2013;Kelly et al., 2019).
TABLE 4 Percent occurrence (Log10 + 1 averages) of COI sequence numbers for 8 taxonomic groups of prey organisms for samples of 7 mesopelagic fish species collected during from 2018 (A) and 2019 (B).
COI metabarcoding can also reveal unexpected prey species in the diets of mesopelagic predators. For example, S. elongatus is known to consume fish and larger invertebrates: two euphausiids, Euphausia krohnii and M. norvegica, were found to be important prey, as well as the squid, Sthenoteuthis pteropus, which was not previously reported. The finding of the fish Benthosema glaciale in the gut contents of C. sloani is consistent with descriptions of the predator as a piscivore, but seven species of euphausiids were also detected by COI. Malacosteus niger has previously been recorded as consuming a variety of crustacean zooplankton, especially copepods (Sutton, 2005); COI sequences identified five fish prey species, excluding the predator species.
Food web dynamics based on molecular data have been examined using ecological networks, where trophic interactions are visualized by pairwise links between species of interest (Clare et al., 2019;Cuff et al., 2022;Dunne et al., 2022), including mesopelagic fish predators and prey (Deagle et al., 2019;Clarke et al., 2020). This study used ecological network analysis of V9 results for the seven predators and 14 prey groups to analyze the relative importance of the many trophic linkages between predators and prey. This view of the food web at the study site shows linkages and interactions that may also apply to other regions and pelagic zones of the Atlantic Ocean. Proportional average numbers of COI sequence numbers (Log10 + 1) for selected taxonomic groups in the gut contents of the 7 predator species collected in 2018 and 2019. No samples of Malacosteus niger (Mn) were analyzed from 2019.
Non-Metric Multi-Dimensional Scaling (NMDS) for COI sequence numbers (Log10 + 1) for 8 prey groups in gut contents of fish collected during both 2018 and 2019. Averages across species include zero values; no samples of Malacosteus niger (Mn) were analyzed from 2019. Samples are identified by abbreviation for the fish species and year of collection.
A fundamental question regarding mesopelagic food web dynamics is whether fish predators exhibit prey selectivity. One approach is to compare biodiversity of the pelagic assemblage with prey composition found in predator gut contents. The relative importance of prey groups in 2018 reported here showed good correspondence with relative abundances of the organisms in the water column, as reported by Govindarajan et al. (2021) based on V9 metabarcoding of MOCNESS samples and environmental DNA (eDNA) collected at the same stations during the HB-1805 cruise. Information on biodiversity of the pelagic assemblage and potential prey groups or species is not available for the 2019 collections during the HB-1907 cruise. In 2018, copepods and ctenophores dominated results of V9 metabarcoding of MOCNESS samples (Govindarajan et al., 2021); these groups also showed high V9 sequence numbers for gut contents of all fish species that we sampled. Hydromedusae and siphonophores were also highly abundant in V9 metabarcoding of eDNA samples (Govindarajan et al., 2021), but these groups were not detected in gut contents.
The selection of target gene regions is an important consideration for metabarcoding: V9 can reliably classify metazoan organisms by taxonomic group, with levels varying among different phyla (deVargas et al., 2015). A 313 bp region of mitochondrial COI (Geller et al., 2013;Leray et al., 2013) has been widely used for metabarcoding analyses of marine metazoan diversity (Hirai et al., 2015;Lindsay et al., 2015;Stefanni et al., 2018;Suter et al., 2021). The shorter COI sequences are somewhat less accurate and reliable for species identification than the COI barcode region of ~675 bp (Bucklin et al., 2011). The eight prey groups analyzed for COI metabarcoding in this study are a subset of the 14 prey groups detected by V9; the missing groups were not detected by COI (i.e., zero sequences for all fish analyzed). The lack of detection across all taxonomic groups of Metazoa was expected,
TABLE 5 Prey counts for selected taxonomic groups based on morphological (microscopic) analysis of gut contents for specimens of the target fish species collected during cruises in 2018 (A) and 2019 (B).
based on previous studies of zooplankton diversity using this COI gene region (Djurhuus et al., 2018;Questel et al., 2021). The analysis of prey using COI metabarcoding must be interpreted cautiously for all metazoan groups, especially including non-crustaceans (Bucklin et al., 2016(Bucklin et al., , 2021a)).
Metabarcoding of any gene region relies on PCR amplification, which can introduce errors, including mismatches between the PCR primer and the target sequence (Taberlet et al., 2012). Also noteworthy is the possibility of DNA degradation, either in ocean waters or in stomach contents of predators. The accuracy and reliability of classification and identification of groups or species based on metabarcoding is also dependent upon the taxonomic completeness and geographic coverage of the reference sequence databases, which must include sequences for specimens identified by morphological taxonomic experts (Bucklin et al., 2011(Bucklin et al., , 2021b)). The recent release and regular updating of the MetaZooGene Database and Atlas (MZGdb) provides increasing confidence in the accurate detection and identification of species and groups across the phylogenetic assemblage of marine organisms (Bucklin et al., 2021b).
Metabarcoding provides a powerful means of tracing trophic pathways and discovering sources of productivity in the deep ocean. Combined use of metabarcoding and microscopy will continue to provide more rigorous and complete reconstruction of mesopelagic food webs than either approach can provide in isolation. The results reported here from integrative analysis using metabarcoding and microscopic counts of gut contents of mesopelagic fishes allowed analysis of primary prey groups and important prey species that sustain the fish species studied. Detailed molecular and morphological taxonomic analysis of net samples from the same depth zones can be used to address questions of prey selectivity.
Future field collections using depth-stratified sampling and carefully-timed deployments to capture both predators and prey with respect to DVM behaviors will be required to determine the extent of selective predation. Environmental differences at the collection sites in 2018 and 2019, including weather and other environmental conditions, may have impacted prey composition of the fish predators. This study was designed to minimize bias by sampling as close to the same time period each year, at the same location defined by geography and bathymetry, using the same methodological approaches, during similar oceanographic and hydrographic conditions. In light of possible ecological, environmental, and oceanographic differences between the 2018 and 2019 sampling events, the results of this study cannot be used to generalize to other seasons or ocean regions. Additional field sampling at regular intervals, with associated environmental and hydrographic data collection, will provide the oceanographic context for prey variability among fish predators, which may result from large-scale, inter-annual or seasonal variation, and/or small-scale, short-term patchiness in pelagic ecosystems of the NW Atlantic Ocean.
The taxonomic completeness and geographic coverage of reference sequence databases are primary limiting factors for detection and identification of taxonomic groups and species based on metabarcoding (Bucklin et al., 2016). Another essential resource is the availability of morphological taxonomic experts for the many metazoan groups comprising the pelagic assemblage (Bucklin et al., 2021a). Several initiatives are currently ongoing that are working toward this shared goal (Duarte et al., 2020;Leray et al., 2022;Bucklin et al., 2021b;Huggett et al., 2022;Govindarajan et al., 2023).
Integrative morphological (microscopic) and molecular (metabarcoding) analysis of gut contents of seven species of mesopelagic fishes was carried out for samples collected during Summer 2018 and 2019 in the NW Atlantic Slope Water. Results are summarized for each predator species, including prey groups detected by V9 and COI, prey species identified by COI, and prey counts based on morphological analysis of gut contents. Patterns of variation for both V9 and COI revealed similarities in diet among the predators, with crustaceans, including primarily copepods and/ or eucarids, with smaller numbers of amphipods and ostracods, dominating in all species except M. niger. V9 allowed detection of gelatinous organisms, including hydrozoans, salps, and ctenophores found in the gut contents of A. aculeatus, C. sloani and S. beanii. The detection and identification of prey species based on COI provided a detailed view of trophic relationships and food web linkages for mesopelagic fish. Parallel analysis of morphological counts in gut contents of fish from the same samples provided additional insights into prey composition. Future research is needed to examine underlying causes of observed differences in prey composition of mesopelagic fish predators, including impacts of vertical migration behavior and whether differences reflect prey choice and availability.
Frontiers in Marine Science frontiersin.org
Bucklin et al. 10.3389/fmars.2024.1411996 Frontiers in Marine Science frontiersin.org