When we initiated this work in late 2023, just three complete ctenophore 28S sequences were in GenBank (Sayers et al., 2023). We used previously published SRA datasets deposited in the National Center for Biotechnology Information (NCBI) to assemble the ribosomal genes (28S) of 50 ctenophores representative of the taxonomic diversity of the phylum as follows ( Supplementary Table S1 ). The first 8 million reads in the dataset were separated and then assembled using the “Map to Reference’’ function and built-in mapper of Geneious Prime 2023 (https://www.geneious.com), with three iterations using an exemplar 28S ctenophore sequence as the seed reference. In contrast to the case of ctenophores, 28S has been used extensively in phylogenetic systematics analyses of medusozoans (Collins et al., 2006; Cartwright et al., 2008; Bayha et al., 2010; Bentlage et al., 2010; Maronna et al., 2016; Miranda et al., 2016; Bentlage et al., 2018). Therefore, to create a local reference database for 28S sequences of Medusozoa, we conducted a comprehensive GenBank search of all available 28S sequences for the clade. The GenBank query used the terms “(large subunit ribosomal rna) AND Hydrozoa[Organism])) NOT Mitochondrial)” and was repeated replacing “Hydrozoa’’ with “Scyphozoa”, “Cubozoa”, and “Staurozoa’’ respectively. The resulting 2,208 sequences were downloaded as a FASTA file.

The Anth-28S-eDNA primer set was tested in silico on the ctenophore and medusozoan 28S reference datasets to evaluate primer complementarity. The bioinformatic process of determining complementarity between the primers and the 28S sequences was conducted separately on the ctenophore and medusozoan datasets but consisted of the same workflow. We used cutadapt (Martin, 2011) to extract all ctenophore and medusozoan sequences that matched the Anth-28S-eDNA primers with two or fewer mismatches and contained at least 10 overlapping bases. Two rounds of cutadapt were conducted on the 28S alignments. First, the cutadapt command contained the flag “action=retained” in order to conserve the sequences matching the primers, including the primer-binding region, in the resulting sequence list. The second round of cutadapt trimmed the primers from the resulting sequences. All sequences were aligned in Geneious Prime (version 2023.2.1) using MAFFT (v7.490) (Katoh et al., 2002). Default parameters used for MAFFT multiple sequence alignment consisted of the Auto algorithm, a scoring matrix of 200PAM/k=2, a gap open penalty of 1.53, and an offset value of 0.123. Forward and reverse Anth-28S-eDNA primers were added to the alignment using the “Add primers to sequence” function with a modified version of Primer3 (version 2.3.7), allowing for two mismatches in the binding region. The resulting annotation table containing base mismatch information was used to compile a summary table with R (version 4.3.3), RStudio (version 2023.12.1 + 402), and tidyverse (Wickham et al., 2019) to visualize the number and location of mismatches between the primers and each 28S sequence. For the ctenophore and medusozoan reference sequences that matched the Anth-28S-eDNA primer set, taxonomic family identifications were added from GenBank.

To test the ability of the 28S barcode amplified using the primers to distinguish between taxa, separate matrices of sequence similarities were calculated between aligned ctenophore and medusozoan sequences using Geneious Prime. Similarities between ctenophore and medusozoan sequences at the family and genus levels were summarized and analyzed in R using tidyverse packages.

Water samples were collected for eDNA analysis during a research expedition on the R/V Point Sur (owned by the University of Southern Mississippi and operated by the Louisiana Universities Marine Consortium) to survey mesophotic banks (30 to 200 meters deep) and deep-sea coral habitats (> 200 meters deep) in waters off the Gulf coast of the southern USA in August 2021. Briefly, eDNA samples were collected using Niskin bottles at the seafloor during remotely-operated vehicle (ROV) dives and in the water column during conductivity, temperature, depth (CTD) casts. Water samples were taken in at least duplicate, if not triplicate or quadruplicate, during each sampling event. Sample collection, eDNA sample processing, and metabarcoding library preparation are described in detail in McCartin et al. (2024) (Supplementary Table S4). Here, we will include information that is critical to interpret the results in the context of the analysis for ctenophores and medusozoans.

The primer set matched especially well to the alignment of ctenophore 28S sequences, with 79.1% of the sequences perfectly matching the forward primer and 89.6% perfectly matching the reverse primer ( Table 1 ). All of the ctenophore sequences matched the primers with two or fewer mismatches. The medusozoan alignment was less complementary, with 10.4% of the sequences perfectly matching the forward primer and 81.1% perfectly matching the reverse primer. Out of all the medusozoan sequences analyzed, 81.7% matched the primers with two or fewer mismatches. Twenty-six of 33 ctenophore families were represented among the reference sequences with two or fewer mismatches to the primer set. Within Medusozoa, six of eight cubozoan families, five of six staurozoan families, 19 of 23 scyphozoan families and 67 of 133 hydrozoan families were represented by reference sequences with two or fewer mismatches to the Anth-28S-eDNA primer set ( Table 2 ).

For both the ctenophore and medusozoan reference datasets, no two families had the exact same 28S barcode ( Supplementary Table S2 ). There was a maximum percent identity of 99.5% and median of 82.8% for ctenophore families and a maximum percent identity of 98.4% and median of 58.6% for medusozoan families, demonstrating that this barcode can be used to distinguish families within these taxa.

After primer trimming, quality filtering and trimming reads, denoising sequencing errors, merging, and chimera removal, there was an average of 12,429 reads per sample with a standard deviation of 8,866. In total, 1,188 ASVs were recovered by the Anth-28S-eDNA primer set and 766 ASVs remained after curation with LULU. Of these 766, 67 ASVs were classified as Ctenophora, and 190 were classified as Medusozoa. Out of the 67 ctenophore ASVs, 58 were classified to at least the order level, 56 were further classified to the family level, 47 were further classified to the genus level, and 3 to the species level. Out of the 190 medusozoan ASVs, 187 were classified to at least the order level, 136 were further classified to the family level, 113 were further classified to the genus level, and 8 to the species level. Within Medusozoa, only Hydrozoa and Scyphozoa were represented. There were 33 families represented, comprising 11 ctenophore families (out of 25 total families) and 24 medusozoan families, specifically from Hydrozoa and Scyphozoa (out of 170 total families).

ROV samples taken at the seafloor and CTD samples taken within 5.8 meters above the seafloor were significantly different at both Bright Bank (PERMANOVA, p = 0.006) and Viosca Knoll (PERMANOVA, p = 0.016). ROV samples were considered representative of the benthic day community and CTD samples were considered representative of the pelagic night community. The benthic community analysis included 25 ROV samples: four from Stetson Bank, 12 from Bright Bank, three from Viosca Knoll, and six from Green Canyon. The pelagic community analysis included 22 CTD samples: 12 from Bright Bank and ten from Viosca Knoll.

The composition of ctenophore and medusozoan communities varied noticeably depending on whether presence/absence data was considered or if relative sequencing read abundance was considered. When presence/absence data was used to perform a dbRDA, the environmental variables accounted for 37.7% of the variance. Samples taken from the same site appeared to group together, and samples collected by CTD from the same site clustered more tightly together than those collected by ROV ( Figure 2A ). All the environmental variables considered were significant except for conductivity, with depth explaining the most variation between samples (PERMANOVA: p < 0.001, R² = 0.12), followed by sampling method of either ROV or CTD (PERMANOVA: p < 0.001, R² = 0.10) and altitude (PERMANOVA: p < 0.001, R² = 0.05). When sequencing read abundance relative to the total amount of ctenophores and medusozoans was considered, 27.5% of variance between samples was explained by the environmental variables. The clustering of samples by sample type or site was not evident on the dbRDA plot ( Figure 2B ). Only conductivity (PERMANOVA: p < 0.02, R² = 0.07), altitude (PERMANOVA: p < 0.02, R² = 0.07), and depth (PERMANOVA: p < 0.02, R² = 0.06) were significant, but explained little of the variation in the data. Any underlying patterns in organismal or DNA abundance could be confounded by bias in PCR amplification.

In total, 192 ASVs were detected and identified to at least the family level, representing 11 ctenophore and 24 medusozoan families ( Supplementary Table S3 ). All ctenophore families detected consist of holopelagic species ( Figure 3A ), whereas medusozoan families include species with holopelagic, meroplanktonic, and benthic habits ( Figure 3B ). There were 24 families represented in the benthic day (ROV) samples compared to 30 in pelagic night (CTD) samples, with 19 ctenophore and medusozoan families shared between benthic day and pelagic night communities. In order of highest to lowest sequencing read abundance, shared families consist of the ctenophore families Cydippidae, Haeckeliidae, Mertensiidae, Lampoctenidae, Charistephaneidae, Ocyropsidae and Bathocyroidae, and the medusozoan families Rhopalonematidae, Diphyidae, Cuninidae, Erennidae, Halicreatidae, Agalmatidae, Physophoridae, Zancleidae, Aeginidae, Eudendriidae, Hippopodiidae, and Nausithoidae. Twenty-one of 56 ctenophore ASVs and 29 of 136 medusozoan ASVs identified to the family level were detected in both benthic day (ROV) and pelagic night (CTD) samples ( Supplementary Table S3 ).

For the benthic day (ROV) samples, unique families were all medusozoan, representing the families Atorellidae, Porpitidae, Pandeidae, Campulanariidae, and Clytiidae. The families unique to pelagic night (CTD) samples contained a mix of ctenophores (“Ctenocerosidae”, Beroidae, Vampyroctenidae, and Cestidae), and medusozoans (Proboscidactylidae, Ptilocodidae, Solmarisidae, Tetraplatiidae, Tubulariidae, Apolemiidae, and Geryoniidae).

At the benthos during the day, two medusozoan families, Rhopalonematidae and Diphyidae, and two ctenophore families, Mertensiidae and Cuninidae, were consistently present at all sites. Dominant families, or families with the highest relative sequencing read abundance, were identified at each site for both Medusozoa and Ctenophora. Within Medusozoa, Rhopalonematidae was the dominant family at Bright Bank and Green Canyon, Diphyidae at Stetson Bank, and Agalmatidae at Viosca Knoll. Within Ctenophora, Mertensiidae was the dominant family at Stetson Bank and Green Canyon, both Mertensiidae and Cydippidae were most abundant at Bright Bank, and Haeckeliidae was dominant at Viosca Knoll. The five families that were exclusive to the benthic daytime samples were all meroplanktonic medusozoans with benthic hydroid stages.

In the water column at night (CTD), the ctenophore families Cydippidae, Mertensiidae, Ocyropsidae, Haeckeliidae, and Charistephaneidae, and the medusozoan families Diphyidae, Rhopalonematidae, Agalmatidae, Proboscidactylidae, Solmarisidae, Ptilocodiidae, Aeginidae, and Cuninidae were present at every site. Medusozoan and ctenophore communities were both dominated by two main families. Rhopalonematidae and Diphyidae were the most prominent medusozoan taxa at both Bright Bank and Viosca Knoll. Mertensiidae and Cydippidae prevailed over other ctenophores, with most ASVs identified as Cydippidae sampled from Bright Bank and as Mertensiidae from Viosca Knoll. Three meroplanktonic families were found exclusively in the pelagic night samples, but it is not possible to determine whether polyp or medusa stages were being detected.

Alpha diversity of benthic daytime communities sampled by ROV followed the same general pattern regardless of diversity measure used. In order of lowest to highest average diversity, the pattern tended to be Stetson Bank, Green Canyon, Bright Bank, and Viosca Knoll ( Figure 4A ). Bright Bank had higher variability in alpha diversity than the other sites. This same trend was represented in the pelagic night data, as Bright Bank consistently had lower average alpha diversity than Viosca Knoll ( Figure 4B ). Although the species accumulation curves show that none of the sites were sampled enough to adequately represent their diversity ( Figure 5 ), they reflect the same pattern seen in the alpha diversity data from Figure 4 .

Based on the in silico assessment, nearly all ctenophores and more than 80% of medusozoans in the 28S reference datasets are likely to be amplified by the Anth-28S-eDNA primers, assuming they can be amplified with up to 2 mismatches to the primers. The high percentage of complementary ctenophore sequences suggests that the primer set can capture a broad range of ctenophore diversity from eDNA samples. Only one base pair in the primer set would need to be changed in order to perfectly match the ctenophore reference sequences. The results indicate that a lower percentage of medusozoans can be recovered, which may be due to the occurrence of mismatches towards the 3’ end of the primer sequence that may interfere with primer binding. However, the primer set is still expected to amplify a large proportion of accepted families across medusozoan diversity ( Table 1 ), suggesting that the Anth-28S-eDNA primers will aid in exploring significant medusozoan biodiversity that might go undetected otherwise. Using one primer set to target multiple taxa — Anthozoa, Ctenophora, and Medusozoa — streamlines eDNA lab workflow, such that more sequencing data can be obtained and used to address multiple taxa with the same amount of time and effort.

Having a metabarcoding primer set that readily amplifies ctenophore and medusozoan taxa presents an exciting development in GZ research by augmenting researchers’ ability to obtain sequencing data for these understudied taxa in a non-invasive, non-destructive manner, and the potential to uncover more about their community composition and distribution. For example, within our data we detected an ASV for the hydrozoan Porpita, a colonial neustonic form known as a “blue button”. Rather than being detected near the surface, however, Porpita was detected as present in 4 of 6 samples near the bottom (sampled by ROV) from below 500 m at both of the Green Canyon sampling points. While it is known that Porpita produces free-swimming medusae, little has been known about their depth distribution (Helm, 2021). However, a recent study sampled surface waters off the Pacific coast of Mexico and detected larval Porpita colonies as well as planula larvae, documenting this stage for the first time (Santiago-Valentín et al., 2024). The surface waters sampled by Santiago-Valentín and colleagues (2024) did not contain medusae, and the production of planula larvae has still yet to be observed. Our study raises the possibility that Porpita medusae sink to bottom waters to carry out the sexually reproductive phase of their lifecycles, although alternative explanations could include the sinking of fecal matter for some organism that fed upon adult Porpita or the sinking of a senesced animal tissue or body. Either of these alternative interpretations for our unexpected detection of Porpita DNA at 500 m involve links between deep benthic ecosystems with those at the ocean’s surface.

We expect that the majority of sequences identified with this primer set will represent members of the gelatinous zooplankton community. However, several taxa are restricted to the benthos and, therefore, are not technically zooplankton (e.g., some groups within Leptothecata and Anthoathecata). Thus, these primers will likely enrich a diversity of medusozoans from the water column to the benthos. Indeed, we detected benthic taxa (e.g. Eudendriidae) in our metabarcode profiles. However, the medusa stages of medusozoans tend to be considerably larger in body size and more dispersed than the polyp stages, which should skew our sampling toward medusae and thus the GZ component of Medusozoa. Other studies have suggested that larger animals are easier to detect in eDNA samples (Derycke et al., 2021).

While 28S is a promising marker for metabarcoding analyses to characterize Ctenophora and Medusozoa assemblages, its power to distinguish species is still unknown. While 28S has been applied as a barcode within sponges (e.g., Erpenbeck et al., 2016; Voigt and Wörheide, 2016; Itskovich et al., 2022), it has been little explored as a viable barcode for ctenophores or medusozoans. Within ctenophores, studies aimed at elucidating species showed mixed success for 28S as a barcode, with some species having distinct sequence signatures but other more recently diverged species being invariant in 28S (Johnson et al., 2022; Alamaru et al., 2017). Regarding medusozoans, 28S has been used extensively in phylogenetic studies of different clades (e.g., Collins et al., 2006; Cartwright et al., 2008; Bayha et al., 2010; Bentlage et al., 2010; Maronna et al., 2016; Miranda et al., 2016; Bentlage et al., 2018) and the 28S marker has been sampled from multiple exemplars of the same species in a few systematics studies (Montano et al., 2015; Cunha et al., 2017). We are not aware of any study that has explicitly assessed 28S as a barcode for medusozoans, but these studies mirror the case with ctenophores, with many distinct species within a genus being differentiated by 28S, but also with some more recently diverged species being invariant in 28S.

In our community analyses, family data were used because the majority of the recovered ASVs would have been excluded if we only used those that were identified to the species level (8 of 257 ASVs). Family was the lowest taxonomic rank that would represent the majority of recovered ASVs (188 out of 257) while still providing useful ecological and biological information.

The two dbRDAs demonstrated the power of metabarcoding with the Anth-28S-eDNA primer set to distinguish medusozoan and ctenophore communities by geographic location and ecosystem type (benthic/day-ROV or pelagic/night-CTD). In addition, using this primer set in conjunction with environmental data can help determine which environmental factors may be responsible for driving similarities or differences between those communities.

GZ assemblages were site-specific, although there was some overlap in community members ( Figure 2A ). The PERMANOVA test results also indicate that GZ communities are structured by both depth and proximity to the seafloor; in both dbRDA analyses, depth and altitude were among the top three variables that significantly explained the most data variance. This assertion can be seen in the tighter clustering of pelagic night samples away from benthic day samples. Previous studies have also found that GZ assemblages are spatially differentiated from one another. Zaldua-Mendizabal et al. (2021) found that GZ assemblages were distinctly differentiated across sites as a result of location-specific environmental conditions, with water column depth as one of the main factors determining assemblage differences. It has also been suggested that hydrographic boundaries and water mass dynamics, which are driven by the physical topography of a particular site, distinguish zooplankton assemblages (Haberlin et al., 2019; Doyle et al., 2007; Youngbluth et al., 2008). The influence of location-specific environmental factors and water column depth may explain why we see spatially separate jelly communities in our results. Geographic distance was not incorporated in the dbRDA analysis and we therefore cannot assert that latitude and longitude are influencing our results. However, we note that sites for which we detected more similar communities ( Figure 2A ) have less geographic distance between them ( Figure 1A ). There may therefore be some correlation between geographic distance and GZ community composition, with higher community similarity between sites that are closer together (O’Donnell et al., 2017).

eDNA transport could influence the interpretation of our results across sites and depths, since currents and mixing can transport particles away from their sources. We consider this possibility, but also note that results from a number of studies suggest that eDNA is distinct across depths in the absence of water column mixing (Andruszkiewicz et al., 2019; Allan et al., 2021; Canals et al., 2021; Govindarajan et al., 2021; Monuki et al., 2021; Govindarajan et al., 2023). Some of these studies analyzed eDNA samples collected with Niskin bottles (Canals et al., 2021; Govindarajan et al., 2021, Govindarajan et al., 2023). In August at our sampling locations, a thermocline was present below depths of approximately 20 meters. This is consistent with limited mixing typical of the warm, late summer at our stations ( Supplementary Figure S1 ). Further, we found that the amount of eDNA from corals in these samples decreased rapidly with altitude above the seafloor, suggesting limited vertical transport near the benthos (McCartin et al., 2024). eDNA metabarcoding data also discriminates between biological communities connected by water movement (e.g. Port et al., 2016; Jeunen et al., 2019). Our field sites are distant and bathymetrically isolated from one another and also represent substantial geological features (e.g. Bright Bank). We can reasonably assume that dispersal of eDNA between our sites is not important for interpreting differences in a broad sense among our sites.

Compared to Figure 1A , there were no clear visual patterns in the data shown in Figure 2B , which accounted for relative sequencing read abundance. Samples do not appear to be grouped by site, ecosystem/sampling method/time of day, or any other environmental variable. Fewer variables are significant — only conductivity, altitude, and depth — than with presence/absence data. This more haphazard grouping suggests that the amount of DNA from a species detected in one sample is highly variable when compared to other samples, even those from the same location or ecosystem type. This high degree of variability, when relative abundance is considered, could reflect local patchiness of DNA shed from individual organisms over small distances, possibly related to the mobility of most medusozoans and ctenophores. Preliminary research on jelly eDNA suggests that eDNA release rate is highly inconsistent from individual to individual (Minamoto et al., 2017). Andruszkiewicz Allan et al (2021) also reported high variability in the concentrations of eDNA sampled from tanks of scyphozoans and suggested that this patchy signaling could be explained by scyphozoans shedding genetic material in larger chunks, since their main body is formed by acellular mesoglea and their cells may not slough off like other organisms. This is consistent with our results, since uneven signals could be caused by jellies sporadically shedding larger eDNA fragments.

A diverse set of medusozoan and ctenophore families were readily detectable using these primers, opening the door to understanding more about GZ community composition and how it varies across space and time within marine environments.

In this study, there are some limitations to assertions we can make due to our sampling. The samples we used were taken during the day at close proximity to the bottom by ROV (at four sites) and in the water column at night by CTD rosette (at two sites). Further, the volume of the samples that was filtered differed between the two methods, because the Niskin bottles on the ROV had a smaller capacity than the Niskin bottles on the CTD rosette. This difference could translate to differences in the spatiotemporal resolution of the methods and also influence estimates of alpha diversity. While sample volume can influence measurements of alpha diversity, samples that differ in volume have been shown to elucidate similar patterns in community composition (Govindarajan et al., 2022; Peres and Bracken-Grissom, 2024). Our data show that GZ communities sampled by ROV and CTD significantly differ at the two sites where both types of sampling were conducted. CTD/pelagic night samples are richer in terms of families (30 vs. 24) and ASVs (162 vs. 80) detected than the ROV/benthic day samples ( Supplementary Table S3 ). One might predict that GZ communities are richer at night than they are during the day due to diel vertical migration (as well as horizontal transport), which occurs in both deep and shallow waters (Harvey et al., 2009; Williamson et al., 2011), when many pelagic zooplankton species make their nocturnal journey to shallower depths in order to feed. However, it is not possible to disentangle whether these differences are due to proximity to the bottom, time of day, sampling method, or a combination of these factors. Future studies with sampling both at the bottom and throughout the water column during the same time frame would more appropriately address what factors drive community differences.