*2.1. Nematode Sampling and Identification*

Marine meiofaunal nematodes were collected in March and September 2011 from the intertidal zone in seven sites along the coast of Quintana Roo State, Mexico (Figure 1, Table 1). Four sediment samples were collected from each site using a Falcon corer (10 × 2 cm) [56]. Individuals were extracted by decantation in the field using two sieves (180 and 63 μm mesh) and fixed with DESS solution [57] or Formalin 10% [58]. In the laboratory, individuals were separated one-by-one and picked up under a stereomicroscope (NIKON SMZ-1). Then, nematodes specimens were individually transferred to temporary slides in a drop of MilliQ water covered by a coverslip and observed with an OLYMPUS BX51 microscope at different magnifications (10×, 40×, and 100×). Well-preserved specimens were identified morphologically and photographed with a camera (Canon G11). When was possible, more than one representative for each morphotype was subsequently selected for further molecular analyses. Morphological identification was carried out using available taxonomic keys for marine nematodes and comparison with original descriptions using several morphological parameters: length and maximum body width, size and position of setae, size and position of amphids, cuticle ornamentation, size and shape of spicula, presence of precloacal supplements, type and tail size and de Man's ratios (a, b and c): a = body length/body width, b = body length/esophagus length and c= body length/tail length [59–61].

**Figure 1.** Sampling localities from marine nematodes in the Mexican Caribbean. Site IDs (A, B, C, D, E, F, and G) are explained in Table 1.


**Table 1.** List of sites where samples were collected. Coordinates are express in decimals \*.

\* Lat, latitude; Long, longitude.

#### *2.2. DNA Extraction, PCR Amplification, and DNA Sequencing*

Under the microscope, each selected nematode was removed from the temporary slide using a fine paintbrush, cut in pieces with a scalpel, and preserved in DESS. The DNA was extracted from single individual animals with a HotSHOT technique [62]. DNA was stored at 4 ◦C for further amplification of the COI gene. PCR reactions were performed in a final volume of 12.5 μL containing 6.25 μL of trehalose 10%, 2 μL of ddH20, 1.25 μL of 10X PCR Buffer, 0.625 μL of MgCl2 50 Mm, 0.125 μL of each primer (10 μM), 0.0625 μL of dNTPs (10 mM), 0.06 μL of Platinum® Taq Polymerase (Invitrogen, Carlsbad, CA, USA), and 2.0 μL of DNA template [63]. We used two primer sets according to [64]: LC01490\_t1 and HC02198\_t1 (TGTAAAACGACGGCCAGTGGTCAACAAATCATAAAGATATTGG/ CAGGAAACAGCTATGACTAAACTTCAGGGTGACCAAAAAATCA); LC01490 and HC02198 (GGTCAACAAATCATAAAGATATTGG/ TAAACTTCAGGGTGACCAAAAAATCA). Thermocycler conditions, for both primer sets, were: 1 min at 94 ◦C, 5 cycles of 40 s at 94 ◦C, 40 s at 45 ◦C, and 1 min at 72 ◦C, followed by 35 cycles of 40 s at 94 ◦C, 40 s at 51 ◦C and 1 min at 72 ◦C, and a final extension of 5 min at 72 ◦C. PCR products were visualized on 2% agarose gels stained with ethidium bromide. Amplicons for LC01490\_t1/LC01490\_t1 were bidirectionally sequenced with M13F and M13R primers [65]; ABI 3730 capillary sequencer (ABI, Thermo Fisher Scientific, Carlsbad, CA, USA) using the BigDye© Terminator v.3.1 Cycle Sequencing Kit (Applied Biosystems, Foster, CA, USA). Sequences were obtained at the Canadian Centre for DNA Barcoding (CCDB) following standard protocols for high-volume samples [66].

## *2.3. Data Analysis*

Sequences were assembled and edited with Codon Code Aligner v 3.0.3 software. The Clustal W program was used for the sequence alignment with default parameters. Phred score [67] was used to assess the quality of the sequences applying the following categories: no sequences = failed, mean Phred < 30 = low quality, 30 < mean Phred < 40 = medium quality, and mean Phred > 40 = high quality. Sequences with low quality and cross-contamination were removed for the analysis. COI sequences are available as part of the project FMN (Free-Living Marine Nematodes from Quintana Roo, Mexico) on Barcode of Life Data Systems (BOLD, www.boldsystems.org) [68].

Genetic divergence was calculated in MEGA v6 [69] using the Kimura two-parameter (K2P) distance model [70]. The presence of stop codons and indels was verified to discard any contaminants such as NUMTs (nuclear mitochondrial DNA segments) [71,72]. Basic Local Alignment Search Tool (BLAST) [73] and Identification Request, on GenBank and BOLD, respectively, were used to identify matches to the DNA sequences generated in this study.

#### *2.4. Phylogenetic Tree Reconstruction and DNA Taxonomy*

A Maximum-Likelihood (ML) tree was reconstructed in MEGA v6 software using 1000 bootstrap replications. The tree was reconstructed with sequences generated in this study and others retrieved from GenBank according to length, quality (stop codons presence), position (concerning barcoding region), and taxonomy (Accession numbers are in Supplement Table S1). The best-fitting substitution model was the General Time Reversible model with nonuniform evolutionary rates and invariant sites (GTR+G+I) and as chosen in MEGA. The COI sequence *Bursaphelenchus* sp. (order Tylenchina) was selected as the outgroup [74].

To delimit evolutionary independent entities of marine nematode species and test the resolution of the sequenced COI gene, we applied three methods: (1) Automatic Barcode Gap Discovery (ABGD) [53] with the following parameters: 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] as distance metrics, (2) Barcode Index Number (BIN) system [54], and (3) mPTP by selecting single-locus species delimitation with *p*-value 0.001 (http://mptp.h-its.org/#/tree) [55].
