**3. Results**

### *3.1. Morphological Identification*

A total of 30 morphotypes from 107 individuals were identified. The list of morphotypes and sampling locations is available in Table 2. The order Enoplida was represented by the highest number of morphotypes (11), followed by Desmodorida and Chromadorida (7), Monhysterida (3), Aerolaimida (1) and Plectida (Table 2). Families Desmodoridae and Chromadoridae were the best represented with five genera each and six and five species, respectively. Eighteen specimens, all morphological very similar, belonged to the *Spirinia* genus. However, based on de Man´s ratios (different proportions of bulb and tail sizes), we considered the presence of two distinct morphotypes: *Spirinia* sp.1 and *S.* sp.2 (Figure 2). Such morphological differences are supported by the results of integrative taxonomy (see below).

In Tulum, one individual resembling the species *Pontonema simile* (Figure 3) was identified. This species was initially described by Southern as *Oncholaimus similis* [76] and redescribed by Filipjev [77] as *P. simile*. The *Pontomema* cf. *simile* reported in the present study had a shorter size (2100 μm) and a shorter spicule (32 μm) compared to the species described by Filipjev [77].

**Figure 2.** *Spirinia* sp.2. Adult male. (**A**) Head showing circular amphid. (**B**) Posterior part of the body showing tail and spicule. (**C**) Anterior part of the body showing esophagus and esophageal bulb. (**D**) Posterior region showing spicule, gubernaculum, and tail shape.

**Table 2.** Morphotypes identified at the collection sites. Species are listed by nematode Order. The number of individuals and COI sequences analyzed in this study are indicated. Taxonomy classification is based on Nemys database (www.marinespecies.org) [12]. Sites are indicated as below: Cancún = A; Puerto Morelos= B; Cozumel, Playa Azul = C; Cozumel, Punta Sur = D; Tulum = E; Mahahual = F; and Xcalak = G.


Among the sampled localities, Cozumel harbored the highest species richness (12 species), which all are new for this locality. From the other localities, we found fewer species but a greater abundance of individuals. *Enoploides gryphus* (Wieser & Hopper, 1967) is a new record for the Mexican Caribbean while *Haliplectus bickneri* (Chitwood, 1956), *Proplatycoma fleurdelis* (Hope, 1988), *Pontonema* cf. *simile* (Southern, 1914) Filipjev, 1921, *Rhips* sp.1, *Epsilonema* sp.1, *Metachromadora* sp.1, *Spirinia* sp.2, *Odonthophora bermudensis* (Jensen & Gerlach, 1976), *Onyx litorale* (Schulz, 1938), *Epacanthion* sp.1, *Halalaimus* sp.1, and *Metaparoncholaimus* sp.1, are new for Mexico.

**Figure 3.** *Pontonema* cf. *simile* (Southern, 1914) Filipjev, 1921. Collected in Tulum locality. Adult male. (**A**) Habitus. (**B**) Head showing dorsal tooth. (**C**) Tail short and rounded. (**D**) Posterior part of the body showing spicule and tail.

## *3.2. Amplification and Sequencing*

Overall, DNA amplification was successful for 67 individuals. However, sequences with low quality were discarded, and 57 sequences from 20 morphological species were considered (Table 2). Individuals of the species *Actinonema* sp.1, *Epacanthion* sp.1, *Metaparoncholaimus* sp.1, *Odontophora bermudensis*, *Oncholaimus* sp.1, *Prochromadorella* sp.1, and some specimens of chromadorids could not be sequenced (Table 2). All COI sequences produced in this study are new for the molecular databases GenBank and BOLD. According to Phred score [67], all sequences were of high quality (mean Phred > 40), longer than 500 bp (502–667 bp), and without internal stop codons. BLAST match values were between 71% and 85% with Nematoda. The mean genetic divergence was 0.43% and 26.45% for intraspecific and interspecific variation, respectively. In most cases, we observed a clear barcoding gap, although low interspecific variation was detected in a few cases (Figure 4). This study provides the first record of mtCOI sequence for *Monoposthia mirabilis*, *Onyx litorale*, *Enoploides gryphus*, *Pontonema* cf. *simile*, *Proplatycoma fleurdelis*, *Xyala striata*, and *Haliplectus bickneri*, and other higher-ranked taxa (Table 2).

### *3.3. Phylogenetic Analysis and DNA Taxonomy*

Maximum-likelihood analysis using the 57 COI sequences obtained in this work and an additional 70 [47,74,78,79] sequences from GenBank showed congruence between DNA barcode and morphological identification. Sequences were grouped in clusters represented by the morphological identity. The tree was supported in recent clades with bootstrap values >90% (Figure 5). The maximum genetic distance values were observed in *Epsilonema* sp.1 clade (4.02%) and *Enoploides* sp.2 (2.24%). The sequences of *Spirinia* sp.1 and *S.* sp.2, with interspecific divergence values of 7.7%, formed two distinct clades.

**Figure 4.** Relative frequencies of Kimura two-parameter (K2P) distances in the COI sequences of marine nematodes analyzed in this study. Frequencies within species are marked in black, among congeneric species in white, and between species from different genera in grey.

**Figure 5.** Maximum-Likelihood tree reconstructed from COI DNA sequences based on GTR+G+I model. Sequences of marine nematodes from the present study are annotated according to their morphological identification. The 20 haplotypes generated in this study are highlighted in blue color. Numbers in parentheses refer to the number of specimens sequenced. Reference of 70 sequences obtained from GenBank are indicated with an asterisk. Accession numbers are in Supplementary Table S1. Higher taxa (families) are indicated behind the line. Classification is based on Nemys database in www.marinespecies.org [12]. Arely Martínez 2020 ©

To ascertain the number of entities, sequences obtained in this work were analyzed using three different methods: ABGD, mPTP, and BINs (Table 3). ABGD analysis with K2P showed 20 initial and 25 recursive partitions with prior maximal distances (PMD) of 0.0010; 20 initial and 21 recursive partitions with PMDs ranging from 0.0017 to 0.0028; and 20 initial and 20 recursive partitions with PMDs ranging from 0.0046 to 0.0129. With the JC69 model, the same results were observed (see Table 3). The mPTP method based on our phylogenetic tree recovered 20 evolutionary independent entities (*p* = 0.001). BINs split *Enoploides* sp.2 and *Epsilonema* sp.1: BOLD: AAU8181 (n = 1) and BOLD: AAU8183 (n = 1) for the first species, and BOLD: AAU8184 (n = 4); BOLD: AAV1242 (n = 1) for the second species. A total of 22 BINs were recovered. Interestingly, results obtained with integrative taxonomy supported that *Spirinia* sp.1 and *S.* sp.2 are two distinct species.

**Table 3.** Evolutionary independent entities recovered from the COI sequences with the three different delimitation algorithms. Relative gap (X) of 1.1, minimal intraspecific distance (Pmin) of 0.001, maximal intraspecific distance (Pmax) of 0.1, K2P [70], and JC69 (Jukes-Cantor) [75] were selected as distance metrics. I = Initial partition, R= Recursive partition. EiEs = Evolutionary Independent Entities.


#### **4. Discussion**

#### *4.1. Diversity and Distribution of Meiofaunal Marine Nematodes in Mexico*

This study contributes to the current knowledge of the diversity and distribution of meiofaunal marine nematodes in Mexico. The faunistic information in Mexico, and Central America, in general, is traditionally scarce, reflecting the low number of taxonomic experts in the area [80]. Until now, a total of 183 species of marine nematodes are known for Mexico, of which seven have type locality in Mexico [16,22,26,81,82]. In the Mexican Caribbean specifically, there are few ecological and taxonomic studies focusing on marine nematodes [17,18,20] with Isla Mujeres and Banco Chinchorro being the most investigated localities. The Chromadorida order was the most represented in both sites [18,20]. Our results confirm the presence of *Terschellingia longicaudata* (de Man, 1907), *Monoposthia mirabilis* (Schulz, 1932), and *Xyala striata* (Cobb, 1920), all species that have been previously reported for Isla Mujeres [20].

This study adds twelve morphological records to the Mexican list of marine nematodes. Specifically, *Haliplectus bickneri*, *Proplatycoma fleurdelis*, *Pontonema* cf. *simile*, *Rhips* sp.1, *Epsilonema* sp.1, *Metachromadora* sp.1, *Spirinia* sp.2, *Odonthophora bermudensis*, *Onyx litorale*, *Epacanthion* sp.1, *Halalaimus* sp.1, and *Metaparoncholaimus* sp.1 are all new records for the country. Due to specimen immaturity and small sample size, it is necessary to collect additional specimens to obtain proper identification for some of the sampled taxa (e.g., Anticomid, chromadorid, and monhysterid).

The report of Desmodoridae taxon is particularly relevant. This family is globally represented with a high number of species (320 species) [80]. Nevertheless, genetic sequence records for this family are scarce compared with the number of species described [80]. Here, we add 20 COI new sequences from Desmodoridae taxa. However, species identification based only on morphology was particularly tricky due to the lack of well-defined diagnostic traits, weak taxonomic keys reference, and absence of updated databases, including species lists [43,47,80].

We report for the first time meiofaunal nematode species from the Cozumel Island, namely *Chromadorita* sp.1, *Catanema* sp.1, *Epsilonema* sp.1, *Epacanthion* sp.1, *Eurystomina* sp.1, two species of *Enoploides*, *Monoposthia mirabilis*, *Odontophora bermudensis*, *Onyx litorale*, *Proplatycoma fleurdelis*, and *Xyala striata*. Cozumel was represented by the highest number of species, supporting that islands, with a high grade of endemism [83], are particularly relevant areas for the study of biodiversity also for meiofaunal species.

Regarding *P. simile*, specimens found in this work are very similar to the ones described by Southern, 1914 [76], due to comparable buccal cavities with a dorsal tooth, tail shapes, and spicules size. However, we hypothesize that the Mexican record may correspond to a new species for science. First, the original record was found in the intertidal zona of West Ireland, far away from our sample locality. Second, spicule size and total length are appreciably different in individuals from the two populations. Additional investigations comparing specimens from both localities are needed to support our hypothesis.

Individuals of *Monoposthia mirabilis*, *Epsilonema* sp.1, and *Onyx litorale*, were respectively represented by the same COI haplotypes even across different sampling sites, including Cozumel Island and continental sites such as Cancún, Puerto Morelos, Tulum, Mahahual, and Xcalak. This observation suggests that such species have a widespread distribution and a low genetic divergence within and among populations. Therefore, individuals of *M. mirabilis* reported for Isla Mujeres [20] and *Onyx litorale* for Cuba [80], could belong to the same, widely distributed, species.

### *4.2. COI Amplification*

We recommend using HotSHOT technique to extract DNA from marine nematodes [62] because, at least in our study, it provided fairly good results. However, the amplification success rate was difficult (85.07%), confirming previous studies [44,47,79,84]. The problem lies in the high mitochondrial genome diversity among nematodes [85–87] and, consequently, a low success rate using universal primers as well as low availability of specific primers for the Nematoda group [47,88]. In five individuals, the length of COI sequences varied (approximately 150 bp) but these sequences were not found related to any nematode taxa at the rank of genus or family. Specifically, we consistently observed a low success of COI amplification for *Epacanthion* sp.1, *Halalaimus* sp.1, *Prochromadorella* sp.1, *Proplatycoma fleurdelis*, *Odontophora bermudensis*, and *Oncholaimus* sp.1, which highlights the need for new primers specifically designed for these taxa. COI sequences obtained by combining existing protocols developed for zooplankton eggs [62] and by Ivanova et al. for invertebrates [63], allowed us to obtain high-quality DNA barcodes and amplify the complete 'barcode region' (> 500 bp). Obtaining long sequences is particularly important; when COI sequences are amplified in a region without a definite barcoding gap (e.g., JB2/JB3) due to overlapping intra/inter-specific distances, the delineation of species boundaries could be unreliable [44,79,84]. Although not all the sequences could be assigned to a species-level taxon, our results will significantly help the taxonomic identification of nematodes in further studies and when high-throughput sequencing approach on environmental DNA is applied.

#### *4.3. Integrative Taxonomy*

COI sequences generated in this work contribute to increasing the record of public genetic repositories and are references for future studies focused on the patterns of diversity and distribution of marine nematodes. In meiofauna and particularly marine nematodes, taxonomic gaps may be compensated only with the integration of DNA and morphological-based taxonomy [31,46]. The combination of both approaches allows us to disentangle diversity at the species level. For example, the two *Spirinias* species showed different morphological traits. However, we could certainly assign them to different species only through integrative taxonomy. In fact, all specimens of *Sprinia* sp.1 were immature, although preliminary observations allowed to ascribe this species to either *S. parasitifera* or *S*. *septentrionalis*, both previously reported for Mexico [20,22]. *Spirinia* sp.2 was represented by several adult specimens. However, none of the morphometric measurements were sufficient to ascribe them to a known species. *Spirinia* sp.2 could be ascribed to *S. amata* or *S. parasitifera* because of the number of precloacal supplements; however, the larger size of spicules and larger total length size suggest

similarities with *S. inaurita*. We do not discard the possibility that *Spirinia* sp.2 could belong to new species to science.

#### *4.4. Application of Species Delineation Models to Disentangle Diversity*

Our results confirm the advantage of disentangling species and estimating marine nematode diversity using the COI gene. Fast-evolving genes and relatively short genes such as mtCOI are not expected to resolve the phylogenetic relationships in the deeper nodes [39,89]. The use of a multiple-gene approach is needed to clarify the tree topology, thereby uncovering phylogeographic relationships and historical biogeography within the group [42]. However, as previous studies suggested, mtCOI can resolve relationships among closely related species [79] and, as suggested by Derycke et al. [47], we support that COI is a useful biomarker to disentangle the diversity of nematodes to the species level. The design of new primers for different taxonomic groups of marine nematodes and the increase of COI data sequences will advance our knowledge of the diversity in different ecosystems [84] as well as provide a better estimate of global species diversity [44].

Analytical methods ABGD and mPTP supported the presence of the 20 species identified morphologically. The ABGD method relies on user-selected parameters (distance model and a prior limit on intraspecific divergence). In our case, we took care not to select a high prior intraspecific divergence (P) to avoid combining all sequences in one single group [53]. For this reason, a prior value was set after analyzing the interspecific divergence values within our generated sequences. The value of 0.01 for prior intraspecific divergence showed the strongest congruence between groups recovered and species defined as proposed by Puillandre et al. [53]. Moreover, we obtained an equal number of groups after the third partition (P = 0.0046) with both models (K2P and JC9); this supports that the threshold for the barcoding gap in the sequences considered in this work is well-defined.

The BINs method showed a 90% match for the recovered evolutionary independent entities. Only *Epsilonema* sp.1 and *Enoploides* sp.2 were split into four BIN numbers. It is essential to consider that this method initially employs single linkage clustering, coupled with a 2.2% threshold to establish preliminary OTU boundaries [54]. If the threshold is higher than 2.2%, the method tends to separate a higher number of entities and overestimate diversity. *Epsilonema* sp.1 and *Enoploides* sp.2 both showed the highest intraspecific divergence values > 4% (4.02% and 2.24% respectively), although they are identical morphologically. Based on all the criteria, we recognized *Epsilonema* sp.1 and *Enoploides* sp.2 as a single species, concordantly to the results obtained by applying ABGD and mPTP methods. We should carefully consider the BIN numbers as a tool for disentangling nematode species until a more robust genetic database is developed.

#### **5. Conclusions**

Our work contributes towards building a more robust taxonomic and genetic database of meiofaunal nematodes, as well as advancing our knowledge of nematode diversity and distribution. At present, DNA sequencing techniques have opened up new possibilities for taxonomic research in meiofaunal nematodes. However, barcode sequences from marine nematodes are underrepresented in light of the diversity of the phylum [34]. Our study combines classical morphology-based taxonomy with DNA sequences to successfully delimit 20 marine nematode species. Although a new set of primers should be designed for some species, our data support that the COI gene represents an excellent molecular marker to disentangle nematode diversity. Moreover, we validate and support the intraspecific distance value threshold of 5% for nematode species [43,47,84]. Lastly, our study reports new records from an unexplored region and contributes to understanding patterns of diversity and distribution of nematodes in Mexico and worldwide.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/1424-2818/12/3/107/s1, Table S1: Records downloaded from GenBank and using for reconstructing the Maximum-likelihood tree.

**Author Contributions:** Conceptualization, A.M.-A. and F.L.; data curation, A.M.-A., A.D.J.-N., and F.L.; formal analysis, A.M.-A. and F.L.; funding acquisition, A.M.-A.; investigation, A.M.-A., A.D.J.-N., and F.L.; methodology, A.M.-A.; project administration, A.M.-A.; software, F.L.; validation, A.D.J.-N. and F.L.; writing – review & editing, A.M.A., A.D.J.-N., and F.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Consejo Nacional de Ciencia y Tecnología (CONACyT) through the Red MEXBOL Project no. 271108.

**Acknowledgments:** We thank Jose Angel Cohuo, Mario Yescas, and Blanca Prado for assistance during fieldwork; Jessica Landeros and Josué Puc for their help during editing; Yareli Cota for nematodes mounts; Holger Weissenberger for elaborating the map; Ruth Gingold (sweepandmore.com) for proofreading earlier versions of the manuscript; and, finally, the Canadian Center of DNA Barcoding (CCBD) for support during the sequencing process.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
