Hidden diversity of Ctenophora revealed by new mitochondrial COI primers and sequences

Abstract The mitochondrial gene cytochrome‐c‐oxidase subunit 1 (COI) is useful in many taxa for phylogenetics, population genetics, metabarcoding, and rapid species identifications. However, the phylum Ctenophora (comb jellies) has historically been difficult to study due to divergent mitochondrial sequences and the corresponding inability to amplify COI with degenerate and standard COI “barcoding” primers. As a result, there are very few COI sequences available for ctenophores, despite over 200 described species in the phylum. Here, we designed new primers and amplified the COI fragment from members of all major groups of ctenophores, including many undescribed species. Phylogenetic analyses of the resulting COI sequences revealed high diversity within many groups that was not evident from more conserved 18S rDNA sequences, in particular among the Lobata (Ctenophora; Tentaculata; Lobata). The COI phylogenetic results also revealed unexpected community structure within the genus Bolinopsis, suggested new species within the genus Bathocyroe, and supported the ecological and morphological differences of some species such as Lampocteis cruentiventer and similar undescribed lobates (Lampocteis sp. “V” stratified by depth, and “A” differentiated by colour). The newly designed primers reported herein provide important tools to enable researchers to illuminate the diversity of ctenophores worldwide via quick molecular identifications, improve the ability to analyse environmental DNA by improving reference libraries and amplifications, and enable a new breadth of population genetic studies.

Pleurobrachia pileus. However, most deep-living ctenophores remain undescribed because specimens are delicate, difficult to access, and are often damaged during collection (Haddock, 2004). The use of remotely operated vehicles (ROVs) and specialized sampling equipment has expanded our ability to observe and collect ctenophores in the deep-sea during the last 30 years .
Morphological investigations of specimens collected from the deep sea suggest that we have only begun to reveal the remarkable diversity within this phylum. However, there are few taxonomic experts on ctenophores, and morphological identifications often are stymied by damaged specimens, cryptic morphology, and poor preservation in all fixatives. Molecular data can provide relatively quick identifications, especially for taxa like ctenophores, although polymorphic loci and good reference libraries are critical to achieve this goal.
Sequence data from the nuclear 18S ribosomal gene provide a molecular phylogenetic framework for the broad relationships within ctenophores (Podar et al., 2001;Simion et al., 2015), and transcriptomes from a handful of species allow for more in-depth studies of some representative diversity (Simion et al., 2017).
However, the 18S rDNA gene fragment is highly conserved, and the phylogenies often do not effectively discriminate between many species and closely related genera, particularly in groups such as Lobata. For example, some genera, such as Bathocyroe, Eurhamphaea, Deiopea and Kiyohimea, have nearly identical 18S rDNA sequences, showing the limitations of the utility of the 18S fragment with respect to species delineation (Haddock et al., 2021).
The "barcoding" mitochondrial cytochrome-c-oxidase subunit-I (COI) sequence fragment is typically useful for species identification and delimitation because it is (usually) easily amplified with universal primers, variable between species, and easy to align between divergent taxa (DeSalle & Goldstein, 2019). Degenerate PCR primer sets are available to amplify COI from many taxa, spanning bacteria to humans (Folmer et al., 1994;Geller et al., 2013;Leray et al., 2013;Siddall et al., 2009). While the COI barcode locus can be successfully amplified for many taxa, it is often difficult for non-model organisms (Vrijenhoek & Waples, 2012). For example, ctenophores have extremely high rates of mitochondrial evolution, are rich in adenine (A) and thymine (T) residues, and have variable gene order in the mitochondrial genome, even within a genus (Arafat et al., 2018;Kohn et al., 2012;Pett et al., 2011;Schultz et al., 2020). A consequence of high mitochondrial variability is that common primers are often poorly suited to amplify the COI fragment. Unsuccessful amplification of ctenophores by commonly used barcoding primers has a number of important ramifications including: a lack of quick and easy molecular identifications that results in difficulties for revealing diversity, few ctenophore sequences in public databases, and a deficiency of easily amplified phylogenetic and population-genetic markers. The lack of COI primers and the resulting paucity of ctenophore sequences available in public repositories also hampers our understanding of the role of ctenophores in ocean ecology. Despite their seemingly low nutrient content, they are frequent prey for a large range of animals Thiebot & McInnes, 2020;Yeh et al., 2020).
Metabarcoding and eDNA studies are powerful tools to assess community diversity, ecosystem monitoring, and function (Eble et al., 2020), but ctenophores are often excluded because they are either missing from reference libraries or are unable to be amplified because of divergent sequences. Some studies have noted the presence of ctenophores in zooplankton metabarcoding analyses based on 18S rDNA fragments (Günther et al., 2018;López-Escardó et al., 2018;Preston et al., 2020;Schroeder et al., 2020;Yeh et al., 2020).
However, since the fragment of 18S lacks resolution to discriminate between most species of ctenophores, often all members of the entire phylum are lumped together, obscuring diversity (Günther et al., 2018;López-Escardó et al., 2018;Preston et al., 2020;Schroeder et al., 2020;Yeh et al., 2020). In studies that used multiple loci in a metabarcoding framework, a significant proportion of 18S sequences were from ctenophores, yet COI often failed to detect any (Djurhuus et al., 2018;Günther et al., 2018), or the ctenophore sequences had poor taxonomic assignments so the results were not informative (Lacoursière-Roussel et al., 2018;Leray & Knowlton, 2015;Pitz et al., 2020).
For this study, we designed multiple primers to amplify COI from all major clades of ctenophores, including many deep-living, undescribed species. We applied those primers across the phylum and tested taxonomic assignments for various groups. Finally, as a case study, we used our newly generated sequences as a library for an eDNA study along the eastern Pacific coast, and provided specieslevel resolution for taxonomic assignments of ctenophores.

| Sample collection
Ctenophores were primarily collected by remotely operated vehicles (ROVs), blue-water SCUBA diving, and midwater trawls, along with contributions from collaborators around the globe (Table S1). Tissue samples for genomic DNA or RNA were flash frozen in liquid nitrogen, then stored at -80°C. Genomic DNA was isolated from frozen samples using the DNeasy DNA Blood and Tissue kit (Qiagen) or the Monarch Genomic DNA Purification Kit (New England Biolabs) according to the manufacturers' instructions. When possible, multiple individuals within a species from different localities and depths were separately amplified and sequenced.

| Primer design
Initial ctenophore COI primers were designed based on COI sequences from our unpublished ctenophore transcriptome data, published mitochondrial sequences of Mnemiopsis leidyi (GenBank accession nos. NC016117, JF760210), Pleurobrachia pileus (GenBank accession no. JF760211), and from the ctenophores which we could successfully amplify with universal Folmer primers (Folmer et al., 1994; Figure 1). RNA extractions, transcriptome sequencing, and analyses were conducted as previously described (Francis et al., 2013). We designed primers for the mitochondrial COI gene with Primer3 (Untergasser et al., 2012)   F I G U R E 1 (a) Primer binding position, taxa amplified, and PCR reaction conditions for each primer set relative to published Mnemiopsis leidyi sequence, including the "Folmer" region and fragment sequences were trimmed to for phylogenetic analyses (b) primer sequences and position used to amplify and sequence all groups of ctenophores including LCO1490/HCO2198 (Folmer et al., 1994) and meta-COI-F2/meta-COI-R2 (Pett et al., 2011) (not used in this study). Differences from M. leidyi are indicated in bold and coloured according to base. Primers marked with * are more universal

| Amplification and sequencing
Fragments of 18S rDNA were amplified using primers MitchA and MitchB (Medlin et al., 1988

| Statistical methods
Sequence fragments were assembled, edited, and aligned with MUSCLE (Edgar, 2004) within Geneious Prime. The COI alignment was translated with the Mold, Protozoan, and Coelenterate translation mitochondrial code table to ensure that the sequences were aligned in the correct reading frame based on annotated GenBank sequences and that no stop codons were present. Where available, we included published sequences for ctenophores for the COI fragment in phylogenetic analyses ( Figure 2).
We used jModelTest (Posada & Crandall, 1998) within Geneious to select the best evolutionary model. We estimated phylogenies with MrBayes (v3.2.7a, Ronquist et al., 2012) and IQ-Tree 2 with 1000 bootstrap replicates (Minh et al., 2020) for all ctenophores available based on the most appropriate models selected by the AIC (Akaike, 1974). Bayesian phylogenies estimated with MrBayes included multiple runs of 5 × 10 6 generations with a 10% burnin, and six chains, that were sampled and printed every 1000 generations.
Phylogenies were visualized with FigTree (v1.4.4,tree.bio.ed.ac. uk/softw are/figtr ee/). To illuminate within-order diversity, alignments of lobate species for both COI and 18S fragments were analysed separately with the same parameters that had been used for the full alignment, and both Bayesian and maximum likelihood support values were reported on the phylogeny. The lobate-specific phylogenies were rooted with the cydippid Pleurobrachia bachei as the outgroup.
To assess diversity within sequenced ctenophores for COI, we calculated percent general time reversible (GTR) distance within and between molecular operational taxonomic units (MOTUs). We considered a clade a MOTU when it was ≤4% GTR distance. Base composition was calculated with MegaX (v10.0.5, Kumar et al., 2018) and pairwise distances were calculated in Geneious Prime.
We calculated the number of parsimony-informative sites within the Lobata for both COI and 18S with DnaSP (v6.12.01, Rozas et al., 2017). We tested for saturation of observed proportions of transitions and transversions versus GTR distance among all MOTUs with the software DAMBE7 (v7, Xia, 2018). We also examined saturation, transitions, and transversions in five other loci from the published mitochondrial genomes of eight species of ctenophores with DAMBE7 to determine if the COI fragment was an outlier.

| Environmental samples and metabarcoding
In 2012 on a transit from Monterey Bay, CA, USA to the mouth of the Gulf of California, Mexico, Pitz et al. (2020) sampled 15 offshore stations with 100 m-depth vertical net tows which they sequenced for 18S and COI metabarcoding fragments. Primers mlCOIintF (Leray et al., 2013) and jgHCO2198 (Geller et al., 2013) were used to amplify the metabarcoding libraries for a ~300 bp fragment of COI. The authors found 9105 COI sequences from thirteen of the stations that were assigned only to the family level or higher (K. Pitz, personal communication). To illustrate how a more complete reference library can affect metabarcoding sequence assignments, we queried the COI sequences previously designated as ctenophores by the Banzai pipeline from Pitz et al. (2020) with our ctenophore-specific library.
Data for stations where no ctenophore sequences were recovered were not reported. We plotted original family and new species assignments based on the presence or absence of reads in each sample in Rstudio (Rstudio Team, 2015) with ggplot2 (Wickham, 2016).

| RE SULTS
In addition to 15 species with published data, we successfully amplified and sequenced 174 individuals of 72 distinct molecular operational taxonomic units (MOTUs). We included published species of Coeloplana, Pleurobrachia globosa, Beroe ovata, and Beroe sp. a and sp. n that we either did not have access to or were undescribed species (accession numbers in Figure 2). We defined a MOTU as a sequence or group of sequences that were ≤4% GTR divergent.
Phylogenetic analyses of the resulting COI sequences revealed high diversity and a number of cryptic species within all groups, but especially within Lobata, which had not been evident from previous 18S rDNA sequence analyses. COI and 18S sequence fragments were deposited in GenBank with accession numbers (MW647008-75, MW735695-832, MW804236-67, MZ391814) (Table S1). Taxonoptimized primer sequences and positions are indicated in Figure 1.    for one species also worked for other closely related taxa, such as Bfor259F, which was specific to B. forskalii but also amplified most

| Primer design and combinations
Beroids. For several species, multiple combinations of primer pairs successfully amplified the COI fragment (Figures 1 and 2).

| Phylogenetics and species delimitation
Combined with published data, the COI alignment included 179 sequences representing 77 molecular operational taxonomic units were ~50% A, 21% T, 16% C, and 13% G, and the vast majority of changes were at the third codon position.
In order to examine saturation between and within taxa, we plotted the proportion of transitions and transversions versus GTR distance among all sequences (Figure 3a). For ctenophores sequenced, we found that the proportion of transversions was greater than that of transitions, which is in contrast to the general rule (for vertebrate and invertebrate mitochondrial DNA), where transitions are often two-fold more common than transversions (Xia et al., 1996).
Transitions were saturated for species that were more than ~30% distinct and transversions were saturated for species that were more than ~25% distinct. We saw the same pattern for the other mitochondrial loci that we could align without the presence of stop codons, including COII, COIII, CytB, ND4 and ND5 ( Figure S1). In these genes as well, transversions were more common than transitions, and there were high levels of saturation.
The high levels of saturation between distantly related species ( Figure 3a) were reflected in the poor resolution of deeper level relationships amongst taxa in the Bayesian and maximum likelihood phylogeny trees. Phylogenetic relationships for COI among all the ctenophores were therefore portrayed in a tree without support values ( Figure 2). In order to illustrate levels of diversity within the phylum, we also plotted within and between pairwise GTR distances for all individuals sequenced (Figure 3b). In comparison to other phyla, distances between specimen sequences were relatively large, ranging between from 0%-43%, and showed four peaks. The first peak represented within-species variability, around 0%-4% GTR, although we had few within-species samples for this estimate. The second small peak centred at 10% GTR distance, was represented by closely-related species complexes. The two largest peaks were represented by the majority of our sequencing efforts and occurred at ~17% (within group; pink, blue, or green) and 30% GTR distance among different orders (grey) of ctenophores (Figure 3b).
The 18S alignment of the Lobata included 82 sequences for 29 MOTUs but included only 47 parsimony-informative sites, was in sharp contrast to the high COI variability. We trimmed the fragment to 1780 bp and had a TrN+i+Γ model determined with jModel Test for phylogenetic analyses. Phylogenetic results showed little differentiation within and even between many genera (Figure 4a).
Posterior probabilities and bootstrap values in the 18S tree were low for most nodes between species. Therefore, most phylogenetic relationships in the Lobata were unresolvable with 18S alone (Figure 4a).
Distinguishing taxa using the COI fragment limited to the Lobata was more successful, since the marker is more variable than 18S. The COI fragment had 199 parsimony-informative sites between 27 MOTUs ( Figure 4b). Although deeper level relationships amongst more distantly related taxa were still not well supported and some saturation was evident (Figure 3a), analyses revealed high support for several cryptic species complexes and new MOTUs. Conversely, other species such as Ocyropsis maculata maculata and Beroe cucumis were revealed to be truly cosmopolitan with worldwide or at least oceanbasin-wide distributions (Figure 4b). were only able to assign taxonomy to the family level using the "Banzai pipeline". In order to illustrate the utility and importance of a good reference library, we used our newly assembled ctenophore sequence library and assigned ctenophore sequences from the five nominal families that Pitz et al. (2020) found ten different species.

| Phylogenetic analyses
Newly-designed primers combined with a few existing primer pairs enabled us to amplify ~700-1000 bp of the mitochondrial COI fragment for 77 MOTUs from all groups of ctenophores ( Figure 2). The COI phylogeny illuminated a great deal of species-level diversity that was not evident for the 18S fragment, within and between groups. However, the levels of divergence and saturation for mitochondrial genes among ctenophores are some of the highest in the Metazoa (Lavrov & Pett, 2016). Despite high levels of mitochondrial variation, within-species diversity of ctenophores rarely exceeded 4% (Figure 3b). One example of high within-species variation was B. infundibulum, which spanned a huge geographic range from the Northern Pacific to the Northern Atlantic Ocean with many geographic barriers within its range (Johnson et al., In preparation).
Saturation within the phylum, especially at third codon positions, Curiously, we found a higher proportion of transversions than transitions for the COI fragment among ctenophores (Figure 3a). For protein-coding loci, transitions are often more common than transversions because they usually result in synonymous mutations and involve the exchange of bases of similar shapes (Xia et al., 1996). It is plausible that the proportion of transitions within our data were overly saturated and additional mutations were obscured. The mitochondrial genomes of ctenophores are enriched with A/T residues, which also could contribute to a greater incidence of transversions (purines⇄pyrimidines) than transitions. Plots of other mitochondrial loci also revealed the same pattern, where transitions saturated quickly, often at ~20% GTR distance, and there were higher proportions of transversions ( Figure S1). Despite high levels of saturation, closely related species and groups were well-supported, especially within the Lobata (Figure 4b). Although saturation among divergent taxa resulted in poor phylogenetic resolution, and it was difficult to make strong conclusions based on the COI fragment alone, COI sequences were easily attributable to MOTUs, and proved useful for species delimitation.

| The Lobata
When the COI phylogeny was limited to the lobate ctenophores, which form a monophyletic clade (including Thalassocalyce and Cestids; Podar et al., 2001), phylogenetic resolution was improved.
The COI fragment revealed that many genera were polyphyletic although many relationships had low support so we were unable to  (Table S1)

| Morphological and molecular concordance
Mitochondrial sequencing reinforced the designation of several undescribed species that returned one 18S MOTU, but differed morphologically or ecologically (Figure 4). Although we hesitate to make conclusions about deeper-level phylogenetic relationships because of high degrees of saturation among more distantly related taxa, we found several MOTUs that confirmed morphological distinctions. For example, for species that were unresolved with the  (Figures 4a,b and 6d-f). We could not definitively assign any Bathocyroe MOTUs to any of the three described species, B. fosteri  from the Gulf of Mexico, B. paragaster (Ralph & Kaberry, 1950)

from the South Western
Pacific, and B. longigula (Horita et al., 2011) from shallow waters in Japan. We collected Bathocyroe aff. fosteri "A" (Figure 6d) from the central California coast, which matched the description of B. fosteri.
Bathocyroe aff. longigula (Figure 6f), also from the central California coast, resembled B. longigula somewhat, in having an elongated stomodaeum, but it lacked any coloured spots along the meridional canals and the gut appeared wider than that of B. longigula. In addition, the B. aff. longigula was mainly found from ~2000 to 3000 m depth and B. longigula from Japanese waters was collected at the surface (although deep species are known to be occasionally upwelled to the surface in that region). Off the central California coast, B. aff.
longigula ( Figure 6f) and B. aff. fosteri "A" (Figure 6d) had overlapping depth ranges and distributions, but had small distinctions in gut shape and were the most divergent from one another molecularly (15% GTR). The third MOTU from the central California coast, Bathocyroe aff. fosteri "B" (Figure 6e), had a distinct gut shape and colour from all the other MOTUs, was relatively rare, and was only collected below 3000 m depth. Unfortunately, specimens from type localities of described species were unavailable for sequencing so it is unclear if any of the MOTUs we sequenced match those already described or new to science. However, with our new primers researchers worldwide can now amplify species of Bathocyroe and help to confirm species identifications.
Previous work posited that the large ctenophore Deiopea kaloktenota (Chun, 1880) might be a juvenile form of Kiyohimea usagi, in spite of the morphological differences between the two species (Matsumoto & Robison, 1992). By sequencing the COI fragment we resolved the relationships between Deiopea and Kiyohimea where rather than just one clade, we found three clades from central coast California, including two distinct Deiopea sp. and K. usagi. A single specimen collected offshore of the Hawaiian Islands also was the same MOTU as one of the Deiopea from central coast California, showing this species has high dispersal capabilities and a broad oceanic distribution (Figure 4b). Specimens of D. kaloktenota, which were described from the Mediterranean, and K. aurita (Komai & Tokioka, 1940) from Japan were unavailable for sequencing, so it is unclear again whether one or both species of Deiopea from Monterey Bay are undescribed.

| Population subdivision of Ocyropsis
The high levels of within-species diversity of ctenophores make the COI fragment, and possibly mitochondrial sequencing in general, a good marker for population genetics. Mitochondrial sequencing revealed many morphologically cryptic but genetically distinct species. Within Lobata, phylogenetic relationships were better supported and many cryptic species complexes were revealed for those with both sympatric distributions and allopatric isolation ( Figure 4b).
The genus Ocyropsis was a good example. Currently, Ocyropsis

| eDNA and metabarcoding
In order to explore biodiversity, ecosystem function, and population genetics, a high-resolution marker with a good reference library is  Figure 5).
When net tow data were queried against sequences generated by the new ctenophore primers, we were able to assign thousands of sequences to 10 distinct species of ctenophores, illuminating how diversity changed over two important biogeographic barriers, Point Conception and Punta Eugenia ( Figure 5).
In conclusion, the primers designed herein will enable researchers worldwide to amplify and sequence ctenophores for the mitochondrial COI locus for quick identification, as a population genetic marker, and for metabarcoding studies. With the publication of our sequences amplified from our new primers, researchers can increase our understanding of marine biodiversity. As we continue to sequence a broader diversity of ctenophores, new sequence information will help with species identification and descriptions, and along with genomic data will also provide a better understanding of relationships within Ctenophora.

ACK N OWLED G EM ENTS
We gratefully acknowledge the crew and pilots of the R/V Western Flyer and ROV Doc Ricketts, R/V Rachel Carson and ROV Ventana for their help with collecting these delicate specimens and assistance with blue-water diving. We would like to thank Warren Francis for generating many of our transcriptome data and for providing helpful comments on the manuscript. April Woods and Maile Prickett were instrumental in PCR and sequencing help, and we thank Kyra Schlining and Susan Von Thun for support at sea. We would also like to thank Jacob Winnikoff who created some of the silhouettes for Figure 2, Genome Legacy Centre to increase our sampling breadth. We also graciously acknowledge Katie Pitz, who provided ctenophore sequences from her metabarcoding library. We also benefitted from the constructive input of three anonymous reviewers who significantly improved our manuscript. This work was funded by the David and Lucile Packard Foundation, the National Science Foundation (NSF-DEB 1542679 to SHDH) and the National Institutes of Health (R01-GM087198).

O PE N R E S E A RCH BA D G E S
This article has earned an Open Data badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available in GenBank with accession numbers (MW647008-75, MW735695-832, MW804236-67, MZ391814).

DATA AVA I L A B I L I T Y S TAT E M E N T
Untrimmed fragments of COI and 18S sequence fragments have been deposited in GenBank with accession numbers (MW647008-75, MW735695-832, MW804236-67, MZ391814) (Table S1). Primer sequences are indicated in Figure 1.