Combined eDNA and Acoustic Analysis Reflects Diel Vertical Migration of Mixed Consortia in the Gulf of Mexico

Oceanic diel vertical migration (DVM) constitutes the daily movement of various mesopelagic organisms migrating vertically from depth to feed in shallower waters and return to deeper water during the day. Accurate classification of taxa that participate in DVM remains non-trivial, and there can be discrepancies between methods. DEEPEND consortium (www.deependconsortium.org) scientists have been characterizing the diversity and trophic structure of pelagic communities in the northern Gulf of Mexico (nGoM). Profiling has included scientific echosounders to provide accurate and quantitative estimates of organismal density and timing as well as quantitative net sampling of micronekton. The use of environmental DNA (eDNA) can detect uncultured microbial taxa and the remnants that larger organisms leave behind in the environment. eDNA offers the potential to increase understanding of the DVM and the organisms that participate. Here we used real-time shipboard echosounder data to direct the sampling of eDNA in seawater at various time-points during the ascending and descending DVM. This approach allowed the observation of shifts in eDNA profiles concurrent with the movement of organisms in the DVM as measured by acoustic sensors. Seawater eDNA was sequenced using a high-throughput metabarcoding approach. Additionally, fine-scale acoustic data using an autonomous multifrequency echosounder was collected simultaneously with the eDNA samples and changes in organism density in the water column were compared with changes in eDNA profiles. Our results show distinct shifts in eukaryotic taxa such as copepods, cnidarians, and tunicates, over short timeframes during the DVM. These shifts in eDNA track changes in the depth of sound scattering layers (SSLs) of organisms and the density of organisms around the CTD during eDNA sampling. Dominant taxa in eDNA samples were mostly smaller organisms that may be below the size limit for acoustic detection, while taxa such as teleost fish were much less abundant in eDNA data compared to acoustic data. Overall, these data suggest that eDNA, may be a powerful new tool for understanding the dynamics and composition of the DVM, yet challenges remain to reconcile differences among sampling methodologies.


INTRODUCTION
Diel vertical migration (DVM) is a well-recognized phenomenon among the world's oceans where large numbers of fish, decapods, cephalopods, and many other species migrate vertically through the water column during crepuscular periods (Hays, 2003;Sutton et al., 2017). Generally, the animals participating in this migration will transit from deep ocean habitats (500-1000 meters) up into the epipelagic depths (0-200 meters) at night and return to the depths at dawn. Many micronekton (actively swimming organisms ∼2-20 cm length; Kloser et al., 2009) perform this daily migration in search of prey, but these migrating organisms are also important food resources for larger predators in the pelagic ecosystem. This striking vertical movement of organisms across different pelagic depths provides opportunities to study unique ecological interactions in the pelagic environment including predator-prey interactions and energy transfer between the surface and deep ocean.
Often, the DVM process is measured using echosounders to detect sound scattering layers (SSLs) in the water column. The SSLs are ubiquitous features throughout the world's oceans and are mainly comprised of vertically migrating taxa that serve as primary prey resources for larger predators, like marine mammals (Proud et al., 2017(Proud et al., , 2019. The SSLs are dominant biological features of these vast ecosystems and contribute importantly to the global process of carbon transport and sequestration (Irigoien et al., 2014). The depths at which SSLs occur are dynamic and often dependent on the geographic location, depth of the water column and time of day, which gives rise to well-recognized and remarkable DVM patterns (Klevjer et al., 2016).
While acoustic techniques can provide information about the dynamics of the DVM, they have a limited ability to discern the specific composition of migrating layers of organism. This gap in knowledge can be filled with towed net sampling, though a single net sampling method often cannot adequately capture all types of pelagic organisms (Milligan et al., 2018). For example, gelatinous zooplankton are not well represented in MOCHNESS sampling gear and thus may be underestimated with net sampling. One promising potential solution is the use of environmental DNA (eDNA) to detect the organisms present in the environment (Thomsen and Willerslev, 2015;Ruppert et al., 2019). eDNA consists of genetic material sampled from an environmental source rather than directly from a biological source (Thomsen and Willerslev, 2015). This sampling technique can capture all the cells that are present in the environment, which may allow for the detection and identification of larger organisms through shed biological trace remnants (cells, feces, mucus etc.). This technique when coupled with next-generation metabarcoding sequencing techniques can provide a census of the organisms present in the environment (Bucklin et al., 2016). To date, eDNA has been leveraged to profile a wide variety of environments including freshwater lakes and streams, terrestrial soils, and marine seawater and sediment (e.g., Bista et al., 2017;Cowart et al., 2018). Studies that have utilized eDNA methods in the marine environment have shown the ability to detect a variety of multicellular vertebrate and invertebrate taxa (Foote et al., 2012;Thomsen and Willerslev, 2015;Cowart et al., 2018), and complement more traditional faunal survey techniques.
In the current study, we sampled eDNA during the DVM in the northern Gulf of Mexico (nGoM) to assess whether this technique captures the dynamics of this daily migration and the composition of SSLs. Sampling of eDNA was directed by shipboard acoustic sensors to sample seawater at different timepoints during the DVM, and eDNA data was compared to insitu and shipboard acoustic data to assess whether changes in eDNA communities were related to changes observed by the acoustic sensors.

MATERIALS AND METHODS
Sampling of environmental DNA (eDNA) was directed by realtime multifrequency acoustic data to identify the dominant SSLs and estimate the depth of the vertically migrating organisms (micronekton). All observations described in this study were collected from a drifting research vessel (R/V Pt. Sur) during one DEEPEND cruise (DP05) where a ship-board multifrequency echosounder system characterized the depth profiles of the SSLs and directed in-situ sampling of CTD rosette outfitted with water samplers and an acoustic probe. For the cast at night, CTD 83 (Sea-Bird 911plus), was deployed within the mesopelagic zone just above the rising SSL (Table 1). It was held in that position while the SSL rose and migrated past the CTD. Seawater samples were collected throughout the DVM process. During the morning CTD cast, the SSL migrated from the surface to deeper, so CTD 84 was placed below the SSL and then held in position while the SSL moved to the same depth and then below the CTD. Water samples were taken at each time-point, similar to night casts. We also sampled in the epipelagic zone (∼90 m), but these samples were not taken during the DVM. The sampling scheme is detailed in Table 1.

Acoustic Data Collection
Multifrequency acoustic backscatter data was collected with a calibrated pole-mounted echosounder system (Simrad EK60 and EK80). The transducers were mounted in an enclosed housing and suspended 2.5 m below the water surface (see Boswell et al., 2020 for additional details). Backscatter data were collected simultaneously at four frequencies (18,38,70,and 120 kHz). The 18, 38, and 70 kHz echosounders had sufficient power and signal to noise to permit examination of the migrating layers, however, effects of attenuation limited data quality beyond 350 m for the 120 kHz. The dominant acoustic features detected from the shipboard echosounder were examined in more detail with an autonomous echosounder (Simrad EK80 WBAT) mounted onto the bottom of the CTD rosette frame and deployed as an acoustic probe (Kubilius et al., 2015;Kloser et al., 2016). The WBAT was configured with a 70 kHz and 120 kHz transducer and was programmed to alternate in 100 ping sequences between the two transducers. The hardware limitations of the echosounder prevented simultaneous collection at both frequencies. The onaxis gain was measured by suspending a 38.1 mm tungsten carbide sphere (∼6% cobalt binder) 7 m below the CTD while probing through the water column. Calibration was performed using the Simrad Lobe program within the EK80 software (Demer et al., 2015). Sound speed profiles and absorption coefficient were computed from bin-averaged CTD data using the Ocean Toolbox (McDougall and Barker, 2011). Raw acoustic backscatter data collected from the shipboard echosounder were imported and scrutinized in Echoview (v9) and processed following methods described by D' Elia et al. (2016) and Boswell et al. (2020). Briefly, data from the transducer face to 15 m depth were excluded from the analysis to account for nearfield effects and to eliminate surfaceassociated interference (e.g., bubble sweep down). Data were examined for interference from other shipboard sonar systems (intermittent or spike noise), false bottom, and background noise. False bottoms were manually excluded. Background noise was identified and removed following a modified process described by De Robertis and Higginbottom (2007). The measurements of Nautical Area Scattering Coefficient (NASC; m 2 nmi −2 ) were derived from the echo integral in 500-m along-track x 5m vertical bins (MacLennan et al., 2002). NASC is generally understood to be proportional to the abundance of biological scatterers and serves as a comparable index of organism biomass (Hazen et al., 2009;Zwolinski et al., 2010;Fennell and Rose, 2015). Integrated backscatter was further binned into three depth intervals: 15-200 m (epipelagic), 200-600 m (upper mesopelagic) and 600-1000 m (lower mesopelagic). We then compared the maximum NASC (m 2 nmi −2 ) values with depth as a function of time from the 18 and 38 kHz transducer by applying a dB-difference ( dB 18−38 = dB 18 -dB 38 ) operation to highlight two categories of scatterers within the scattering layer: (1) swim bladdered fish (SBF) and (2) non-swim bladdered fish/crustaceans (NSBFC) following D' Elia et al. (2016). Swim bladdered fish targets were classified as having a dB 18−38 range between +3 and +12 dB while the non-swim bladdered fish/crustaceans were classified as having a dB 18−38 range between -14 and 0 dB. Backscatter from these two classes were integrated into 20 m (vertical) by 3 min (horizontal) cells and exported from Echoview for additional analysis.
Acoustic backscatter data from the WBAT were examined in Echoview. Given that the emphasis was to enumerate the individual scatterers within the migrating layer, we focused the analysis on the detected point targets using a single target detection algorithm in Echoview. Sequential targets which satisfied the criteria for the algorithm were tracked and identified as an individual if they had more than 3 targets across 5 consecutive pings. For point target recognition, we used a minimum uncompensated TS threshold of −75 dB and an echo length ranging from a minimum 0.7 ms to a maximum 1.50 ms. The maximum gain compensation was set to 6.0 dB with a maximum phase deviation of 0.6 degrees for both major and minor axes. Point targets within 2 to 75 m of the WBAT were included for use in analysis. During the deployment, the WBAT remained stationary at two depths during the migration of the scattering layer: ∼92 and 320 m. The acoustic data were divided into two periods and were related to the depth of the CTD during the profiling activity. The first period occurred between 16:45 and 20:00 and targeted single targets at ∼320 m depth. Data collected during this period were divided into 10-min time intervals for processing. The second period occurred between 20:00 and 20:24 and corresponded to when the CTD was positioned at ∼90 m. Data collected during this period were divided into 3-min time intervals to capture single targets at higher temporal resolution. Target strength (dB re 1 m 2 ) distributions were derived from each interval and analyzed to examine the temporal changes in the size of detected point targets (Lurton, 2002). We conducted a Kruskal-Wallis Chi-squared test to test for significant changes among depths.

Seawater Collections for eDNA
Seawater samples were collected according to the position of the CTD relative to the SSLs and at specific times using the 12 L Niskin bottle array also mounted on the CTD (Table 1).
Seawater was collected to sample the eDNA for potential changes in composition of prokaryotic and eukaryotic organisms. Collections and filtering have been described by Easson and Lopez (2019). Immediately after the CTD arrived back on the surface, replicate two-liter samples (N = 3) were taken from each Niskin bottle and were filtered through sterile 0.45 υm filter membranes to remove all cells from the seawater. No organisms were directly sampled, and filters were visually inspected to ensure no organisms were visible on the filter. Filter membranes were placed in sterile tubes and frozen until analysis at NSU.

eDNA Sequencing and Statistical Analyses
DNA was extracted using the PowerLyzer PowerSoil DNA Isolation Kit (QIAGEN) following the manufacturer's protocol. Polymerase chain reaction (PCR) was performed following the 18S Illumina Amplicon protocol of the Earth Microbiome Project 1 . Sample preparation and sequencing followed the methods outlined in Easson and Lopez (2019), except that a 300 cycle Illumina MiSeq kit was used to generate paired-end 150 bp amplicons.
Bioinformatics processing was conducted in R using the DADA2 pipeline (Callahan et al., 2016;R Core Team, 2018). Initially, sequences were trimmed to remove ambiguous bases (max N = 0), sequences longer than 150 bp and shorter than 90 bp. The default ASV2 parametric error model was used to calculate sequence error rates. Next, sequences were dereplicated to infer sequence variants, paired-end reads were combined, and chimeras were removed. Once these processes finished, a sequence table of amplicon sequence variants (ASVs) was constructed. ASVs were taxonomically classified using the Silva database (release 128; Quast et al., 2015).
Statistical analysis was done in R (R Core Team, 2018) using the R packages vegan and picante (Kembel et al., 2010;Oksanen et al., 2017). Analysis was conducted separately on CTD 83 (ascending DVM at night) and CTD 84 (descending DVM descending DVM during day). ASV abundance was transformed to relative abundance so each ASV is represented as a proportion in the whole dataset for a specific sample. Initial analysis assessed eDNA community diversity (Inverse Simpson's index) and richness (total number of unique ASVs) and compared these metrics across sampling time-points in each CTD. Beta diversity was then assessed by calculating Bray-Curtis dissimilarity among samples and a permuted multivariate analysis of variance (PERMANOVA) was used to compare beta diversity among sampling time points in each CTD.

Shipboard Acoustic Data
During the day, the depth of the maximum backscatter (NASC, m 2 nmi −2 ) occurred at a depth of ∼440 m (NASC = 58.86 1 earthmicrobiome.org m 2 nmi −2 ) ( Figure 5). The scattering layer began migrating upwards at ∼17:30 and completed the migration at 20:23, with a large amount of backscatter filling the upper 150 m of the water column (CTD 83). The downward migration began around 04:00 with the descent completing at ∼08:30 where the scattering layer settled at ∼460 m (NASC = 55.33 m 2 nmi −2 ) (CTD 84).
Two distinct scattering layers were identified which contributed significantly to the measured water column backscatter. When the CTD was in the upper mesopelagic zone (16:45 to 20:00), the center of mass of the SBF was ∼325 m (± 97.9 m). In contrast, the center of mass for NSBFC were located an average of 407 m (± 114.3 m). When the CTD was raised to the epipelagic zone (20:00 to 20:25), the center of mass of SBF was 155 m (± 43.6 m), while NSBFC were at ∼292 m (± 217 m). Variability in the depth of the NASC maxima between SBF and NSBFC was consistent during both depth profiles.
We observed significant variability when examining the temporal patterns of target strength distributions collected at the two depths from the probe positioned within the ascending DVM. Additionally, the target strength distributions varied significantly in the distribution of scatterers within the mesopelagic layer at 320 m (p < 0.001; χ 2 = 830.8, df = 17) but not within the epipelagic layer (p = 0.629; χ 2 = 5.253, df = 7) (Figure 1).

eDNA Profiles of nGoM Organisms
A total of 3,033,134 sequence reads were recovered for the 40 samples in CTD83 and CTD84. Taxonomic classification using the Silva database (Release 128) identified 10,185 unique ASVs. Since the goal of this study was to assess detection of micronekton that are potentially vertically migrating, taxa that did not fit into this category based on taxonomic classification or had an unknown classification were excluded. Most taxa in the dataset were either unicellular (3,634 taxa) or could only be definitively classified as eukaryotes (6,326 taxa) and were thus excluded from analysis (Supplementary Table S1 An exception to unicellular taxa not being involved in the DVM may be those unicellular taxa that are parasitic to vertically migrating taxa such as in the order Syndiniales, which were abundant in the current dataset (Guillou et al., 2008). Taxonomic classification of eDNA for most ASVs was only to order due to a short amplicon length (∼150 bp), limited taxonomic resolution in the Silva database, and a lack of many deep-sea organisms in taxonomic databases. In many cases, unique ASVs were detected, but these ASVs could only be classified to a higher taxonomic group. Henceforth, when individual taxa are referenced, these are unique ASVs that may be only broadly taxonomically classified.

CTD 83 -Ascending DVM
For the ascending phase of the DVM in the mesopelagic, we observed significant shifts in eDNA diversity and composition over time (Supplementary File 1) with several groups of organisms exhibiting temporal shifts in relative abundance, suggesting variability in the structure of eukaryotic communities Frontiers in Marine Science | www.frontiersin.org FIGURE 1 | Target strength (dB re 1 m 2 ) density distributions during CTD cast 83 derived from the WBAT deployed at 320 m. Each distribution is approximately representative of a 12 min period where the WBAT was deployed. The asterisks indicate intervals where the water was sampled for eDNA. over time. Broad taxonomic groups such as arthropods, tunicates, annelids, cnidarians, ctenophores, and teleost fish were variable over the course of the DVM (Figure 2). Interestingly, different groups of organisms displayed peak relative abundance at different times during the DVM. For example, copepods in the class Maxillopoda appeared most abundant within the earliest sampled time-points of the DVM (at 19:00) whereas other taxa such as tunicates, ctenophores, and cnidarians appeared most abundant in the subsequent sample taken 26 min later. Teleost fish (Class Actinopterygii), which are typically considered to comprise much of the vertically migrating taxa biomass, did not display a high relative abundance in the eDNA dataset though they were consistently present (Figure 2). During the ascending DVM, gelatinous zooplankton (tunicates, ctenophores, and cnidarians) all exhibited the highest relative abundance within the middle time point of the vertically migrating band (at 19:26), however, since these organisms tend to only weakly scatter sound unless occurring in swarms (Ressler, 2002;D'Elia et al., 2016), it may be difficult to discern them in the acoustic dataset. Conversely, there is another group of gelatinous zooplankton, physonect siphonophores, that do scatter significant amounts of acoustic backscatter due to their gas-containing pneumatophore (Warren et al., 2001;Lavery et al., 2007). It is often not possible to acoustically differentiate scattering contributions from siphonophores and some swim-bladdered fish (Proud et al., 2019) which further complicates the analysis. Lastly, Annelids (Phylum Annelida) represented an abundant group in the eDNA dataset (Figure 2). Annelids in the class Clitellata exhibited an increase in relative abundance late in the ascending DVM (Supplementary Table S2).
The eDNA dataset allows us to breakdown each of these broad taxonomic groups to determine how many potential species are driving the observed dynamics (Supplementary Table S2). Previous research has noted variable vertical migration behavior within broader taxonomic groups, and eDNA techniques could help refine our understanding of these behaviors. Within our current dataset, we observed that in most cases, the main drivers of the DVM patterns (Figure 3) were due to only a few taxa in each group. For example, in Arthropoda, we observed 14 total taxa ( Table 2), but only three of these taxa had much impact on the temporal variation of the group (Figure 3). Two of these taxa (ASV18, ASV19) exhibited a peak relative abundance early in the DVM, while the third taxon, ASV20, was most abundant late in the DVM. Each of these three taxa were copepods belonging to the order Calanoida. Similar to the results in Arthropoda, relatively few taxa were observed driving temporal variation in Cnidaria (3 dominant members from the order Siphonophorae), Ctenophora (1 dominant member), Tunicata (10 taxa with 2 main drivers from the family Oikopleuridae), teleost fishes (1 dominant member), and Annelids (1 dominant member).

CTD84 -Descending DVM
In the mesopelagic zone during the descending vertical migration of organisms, we also observed significant shifts in eDNA diversity and composition over short temporal scales (Supplementary Results). Similar to the ascending DVM in CTD 83, several broad groups of organisms were dynamic over the course of our sampling, however, the patterns were somewhat different in CTD 84. Tunicates in class Appendicularia showed a slight increase early in the vertical migration and then a gradual decrease in relative abundance over time (Figure 4). Arthropods in the class Maxillopoda (likely copepods) decreased in relative abundance until the middle of the vertical migration and then gradually increased in relative abundance. Ctenophore and Cnidarian relative abundance both peaked during the middle of the DVM. Teleost fish displayed low relative abundance in this dataset, but they were consistently present with high variability among samples. Annelids in the class Clitellata displayed a higher relative abundance late in the descending DVM (Figure 4), similar to the observed pattern in the ascending DVM (Figure 3).
The taxonomic composition for the descending DVM was somewhat different compared to the ascending DVM (Supplementary Table S3). In Arthropoda (18 taxa total), only ASV20, which appeared later in the DVM in the ascending phase of the DVM (CTD 83), was abundant in the descending DVM (Figure 5). Two other taxa (ASV 24 & ASV 27) in the order Calanoida and one taxon (ASV66) in the order Harpacticoida, were abundant in the descending DVM eDNA dataset. Fewer tunicate taxa were detected in the descending DVM, and we observed an absence of the two most abundant taxa from the ascending DVM. Annelids and teleost fish (Class Actinoptergii) exhibited lower richness compared to CTD83 and were comprised of only two and one unique taxa, respectively. Within the phyla Cnidaria and Ctenophora, we detected 15 and 8 unique taxa, respectively, though both were generally less abundant compared to CTD83, and were dominated by only a few taxa (Figure 5). The majority of Cnidarian taxa belonged to the order Siphonophorae (9 taxa), which includes taxa with gas inclusions (in suborder Physonectae), but taxonomic identification was only possible to order with the current sequence data.

Epipelagic eDNA Communities
Our samples in the epipelagic zone, which occurred during the DVM, represent post-migration samples for both CTD 83 (night) and CTD84 (day). The eDNA in these samples were composed of similar groups of taxa to those in mesopelagic including copepods (Class Maxillopoda), tunicates (Class Appendicularia), Annelids (Class Clitellata), and hydrozoans (Class Hydrozoa; Supplementary Figure S1). Of these six groups, Maxillopoda was the richest, containing 33 unique ASVs. Within this group, some of the same taxa were dominant (e.g., ASV18), but overall, there appeared to be more copepods in the orders Cyclopoida and Harpacticoida compared to samples from the mesopelagic zone. Cnidarians in the class Hydrozoa were less prevalent and less rich in our epipelagic samples, though some similar taxa were detected. The full results for the epipelagic taxa are summarized in Supplementary Table S4.

DISCUSSION
Studying the oceanic pelagic environment is challenging due to remoteness, depth, and the vast expanse of ocean that this zone covers. Given these fundamental challenges, cost-effective, yet comprehensive tools are needed to capture the biodiversity and dynamics of this understudied environment. Remote sensing tools such as echosounders have been utilized to capture the placement and movement (i.e., DVM) of organisms in the pelagic environment (Milligan et al., 2018), but the ability of this technology to resolve the identities of organisms is limited (D'Elia et al., 2016). The current study suggests that eDNA has the potential to fill in some of these gaps in our knowledge through its ability to resolve finer scale taxonomic identities of the organisms in this environment. Our results show that during the DVM, we are able to detect changes in the eDNA community that are concurrent with shifts in the depth of the SSLs through the migration process. The eDNA results show a complex community of arthropods, tunicates, cnidarians, ctenophores, annelids, and fish, which are all known members of the pelagic environment and in some cases are known diel vertical migrators.

eDNA Faunal Assessment
Our eDNA dataset detected many resident fauna of the pelagic environment, and while we had more resolution on organism identity compared to acoustics, we were often not able to Individual amplicon sequence variant (ASV) classifications for CTD 83 and CTD 84 are located in Supplementary Tables S1, S2, respectively.  resolve fine-scale taxonomic identities. With our short 18S rRNA fragment, we were only able to identify ∼31% of our sequences further than a kingdom classification. The presence of many unidentified sequences may represent additional biodiversity in key pelagic groups that we were unable to classify here. Although efforts have been growing, the pelagic environment is understudied relative to many other habitats (Sutton, 2013;de Vargas et al., 2015;Sieracki et al., 2019). The poor classification of many of our sequences may be due to a lack of representation of whole faunal groups or at least an absence of closely related taxa in the Silva database. These instances would lead to low taxonomic resolution in our dataset despite being able to parse individual taxa via our sequenceclustering algorithm. In addition to our classification issues, we also observed the absence or low abundance of some faunal groups that we expected to be more abundant. For example, teleost fish (Class Actinopterygii) were present, but only in low abundance in our dataset. We hypothesized that much of the biomass observed in the acoustic data was attributed to teleost fish, but this dominance is not reflected in the eDNA dataset. We also expected to observe other groups of crustaceans such as decapods in our eDNA dataset (Thomsen and Willerslev, 2015;D'Elia et al., 2016), but surprisingly these groups were absent. Potential explanations for these findings include: (1) low concentrations of eDNA from these organisms present in the environment during the short timeframe we sampled (low rate of eDNA production) or the actual low abundance of these organisms, (2) a lack of resolution in the taxonomic database to identify organisms in these taxonomic groups (organisms are present in dataset but unidentified), (3) organisms like physonect siphonophores were contributing to acoustic backscatter ascribed to teleost fish (Proud et al., 2019), and (4) issues with sequence processing methodology or the presence of inhibitors in certain cells that prevented proper sequencing of taxa in these groups. Previous studies have noted the lack of coverage of many key pelagic taxonomic groups in major databases such as NCBI and Barcode of Life Data System (BoLD) (Kvist, 2013;Cowart et al., 2018), and at present, these limitations likely limit the overall ability of 18S eDNA to finely resolve biodiversity in these hard to study environments. We identified a diverse array of taxa in the pelagic environment that corresponds with previous faunal assemblage characterizations ( Table 2). In our present dataset, sequences matching several groups of tunicates were identified, and taxonomic classification of these sequences indicated that they belonged to the groups Appendicularia (Larvaceans) and Thaliacea (Salps & Pyrosomes). Previous research has shown that members of each of these groups can alter biogeochemical cycling in the pelagic environment, through carbon export to deeper depths (Andersen, 1998). For Appendicularia, the dominant members of our eDNA dataset, this is through their contributions to marine snow via their discarded mucus houses (Barham, 1979), which can be a food source for other pelagic organisms such as copepods (Steinberg, 1995). Representative taxa such as Pyrosoma atlantica have been observed and collected via MOCNESS sampling at depth in the nGoM during DEEPEND cruises , and this species has been identified in DVM elsewhere (Henschke et al., 2019). In addition, our eDNA dataset revealed as many as 17 Cnidarian taxa belonging to the class Hydrozoa. These hydrozoans were classified into three different orders, with Siphonophorae being the most abundant and species rich. These taxa are broadly distributed in all oceans forming crucial trophic links in the deep sea (Robison, 2004;Mapstone, 2014). Several members of this order contain gas-inclusions (suborder Physonectae) and are therefore detectable by our acoustic sensors. Our current eDNA dataset cannot determine what proportion of the community is physonect siphonophores, but previous studies in the GoM have shown this group can be prevalent among siphonophore assemblages (Sanvicente-Añorve et al., 2007). Ctenophores (phylum Ctenophora), while generally less abundant in the eDNA dataset than tunicates or cnidarians, were consistently detected in both the ascending and descending DVM samples. Ctenophores in four orders were present in our dataset along with three taxa that could not be classified past phylum. Some species of ctenophores and siphonophores have been observed to vertically migrate when strong physical stratification of the water column is absent, thus some taxa present in our dataset may represent active migrators over the course of our experimental timeframe (Júnior et al., 2015). Copepods (Class Maxillopoda) appear as the most dominant potential migrator in our eDNA dataset, comprising nearly 20% of eDNA sequences at some time points. Copepods are abundant zooplankton members in many ocean basins, and previous research has observed variable DVM patterns within this group (Bollens and Frost, 1991;Hays et al., 2001). Across our eDNA dataset, we identified at least 20 copepod taxa, of which 13 were classified as Calanoida. The dominance of calanoid copepods is consistent with other research showing that this group is the most successful of all copepods and able to colonize all depths in the water column (e.g., Huys and Boxshall, 1991;Fosshagen et al., 2001).

Signatures of Copepods: Actual Dynamics of the DVM?
The acoustic data in the present study capture the movement of mostly fish and other strong scatterers that form the dominant members of SSLs (> −60 dB). A clear migration signal can be seen for these organisms rising toward the surface displaying the characteristic DVM behavior, which co-occurs with a gradual increase in the size of organisms around the CTD (Figure 5). While our eDNA dataset does not necessarily capture an increase in these large migrators, we do see an increase in the relative abundance of copepod sequences (Figure 6) that precedes the arrival of larger acoustic targets (i.e., micronekton). Previous research has shown that calanoid copepods, the dominant copepod in our eDNA dataset, are often numerically dominant in these habitats and will migrate vertically to avoid predation Frost, 1989, 1991;Proud et al., 2019). While we lack direct evidence for predation avoidance in the current study, the pattern of relative abundance for two calanoid taxa (ASV18 & ASV19) in CTD 83 may suggest that these taxa are ascending to avoid vertically migrating micronekton predators. Further investigation of these dynamics is warranted, including direct sampling of these organisms, to derive any causality for the observed copepod dynamics, since other research has shown copepod vertical migration can be related to body condition (Hays et al., 2001) or environmental cues (Batchelder et al., 2002).
There are a few points of uncertainty in the current study that prevent us from drawing a direct relationship between the eDNA and acoustic datasets. The 18S eDNA dataset estimates a relative abundance of sequences in a volume of water, which can be difficult to relate to absolute abundances captured in the acoustic data. Additionally, we have some uncertainty as to whether our eDNA samples encompass a definitive time frame given the potential for differential residence time (i.e., how quickly DNA degrades in the pelagic environment; Rees et al., 2014;Thomsen and Willerslev, 2015) and how the presence of an organism in the environment relates to the timing of its presence in the eDNA dataset. Previous research has indicated a relatively low residence time for DNA in seawater (Dell'Anno and Corinaldesi, 2004), which would support that our eDNA dataset dynamics are capturing recent shifts in organism abundance. Despite these limitations, our data show the potential for eDNA to be widely used as a tool for monitoring pelagic biodiversity and the dynamics of the DVM. Several aspects of this study could be improved in future studies to increase the impact and resolution of eDNA data, and here we have made a few recommendations. Since the eDNA dataset captured many smaller pelagic organisms, coupling future eDNA samples with net sampling or a higher acoustic frequency that could discern smaller organisms such as copepods, would help calibrate eDNA results. These comparisons could help answer some key questions in eDNA research such as providing estimates of spatial and temporal distributions and more precise links between the presence of organisms in the environment and the presence of their DNA. Using multiple primer sets [e.g., 18S -current study, COI (fish and molluscs), 16S (crustaceans); Bernard et al., 2017;Bracken-Grissom, 2017] in order to minimize primer bias FIGURE 6 | Relative abundance of copepods detected with eDNA from water samples collected at depth (arrows) relative to the deep scattering layer observed at 38 kHz. Two water collection events are displayed (red line, CTD 83 and 84) during the ascending and descending phase of the migration, respectively. The WBAT was deployed on the CTD during the ascending migration phase at 320 m and 93 m depth). The two black vertical regions in the echogram represent periods where the echosounder was not operational.
in PCR reactions and generating longer sequencing amplicons would likely help capture additional diversity in the pelagic environment. Alternatively, a non-PCR metagenomics approach would eliminate all PCR bias but would require greatly increased sequencing depth.

CONCLUSION
In conclusion, we have expanded the application of an established molecular ecology method (eDNA sequencing) to an important global phenomenon (DVM). We demonstrate the identification of a diverse marine community within our eDNA dataset that for certain groups may capture additional biodiversity compared to more traditional sampling methods. Our data indicate that the taxa present before the DVM (outside the SSL) differ from those present within the SSLs and after the DVM. Some dynamics within our dataset appear to mirror changes in the acoustic data suggesting eDNA metabarcoding is capturing real shifts in the composition of vertically migrating organisms, however, additional effort is needed to resolve the source of backscatter data collected at depth and the identity of organisms that are being detected in eDNA profiles during DVM activities.

AUTHOR CONTRIBUTIONS
CE, KB, and JW contributed to the conception and design of the study. JL and KB collected acoustic and seawater samples. CE generated genetic data, organized datasets, performed eDNA analysis, and constructed the initial design for the manuscript. KB, NT, and JW processed acoustic and CTD data for analysis. All authors contributed to the article and approved the submitted version.

FUNDING
This research was made possible by a grant from The Gulf of Mexico Research Initiative to support the consortium research entitled "Deep Pelagic Nekton Dynamics of the Gulf of Mexico" administered by Nova Southeastern University. This was contribution #200 from the Center for Coastal Oceans Research in the Institute of Water and Environment at Florida International University.