Diversity of free-living marine nematodes (Enoplida) from Baja California assessed by integrative taxonomy.

We used morphological and molecular approaches to evaluate the diversity of free-living marine nematodes (order Enoplida) at four coastal sites in the Gulf of California and three on the Pacific coast of Baja California, Mexico. We identified 22 morphological species belonging to six families, of which Thoracostomopsidae and Oncholaimidae were the most diverse. The genus Mesacanthion (Thoracostomopsidae) was the most widespread and diverse. Five allopatric species, genetically and morphologically differentiated, were found in two localities in the Gulf of California (M. sp1 and M. sp2) and three in the Pacific coast (M. sp3, M. sp4 and M. sp5). Overall, we produced 19 and 20 sequences for the 18S and 28S genes, respectively. Neither gene displayed intraspecific polymorphisms, which allowed us to establish that some morphological variation was likely either ontogenetic or due to phenotypic plasticity. Although 18S and 28S phylogenies were topologically congruent (incongruence length difference test, P > 0.05), divergences between species were much higher in the 28S gene. Moreover, this gene possessed a stronger phylogenetic signal to resolve relationships involving Rhabdodemania and Bathylaimus. On the other hand, the close relationship of Pareurystomina (Enchilidiidae) with oncholaimids warrants further study. The 28S sequences (D2D3 domain) may be better suited for DNA barcoding of marine nematodes than those from the 18S rDNA, particularly for differentiating closely related or cryptic species. Finally, our results underline the relevance of adopting an integrative approach encompassing morphological and molecular analyses to improve the assessment of marine nematode diversity and advance their taxonomy.

Introduction

Nematodes are the most dominant and diverse meiofaunal group in marine benthic habitats. Usually, they account for 70–90% of meiobenthic metazoan abundance in marine sediments, where they play fundamental ecological roles (Austen 2004). Estimates of marine nematode diversity are in the millions of species, of which only a small fraction has been described (De Ley et al. 2005; Lambshead and Boucher 2003). Unfortunately, the field of marine nematology remains underdeveloped, perhaps due to labor-intensive techniques, a difficult taxonomy, and a very limited worldwide expertise (De Ley 2000; Bhadury et al. 2006).

Nematode systematics faces the fundamental question of whether morphological characters alone are sufficient to achieve a natural classification of all nematodes and whether these characters are reliable for species identification. The practical limitations and challenges imposed by natural variation in the identification of marine nematodes are many. For instance, taxonomic knowledge and identification keys are restricted to few geographic areas, and a high level of taxonomic expertise is needed to work with them. In addition to human-related limitations, natural variation in the form of phenotypic plasticity and cryptic/sibling speciation also imposes challenges making morphological variation taxonomically equivocal. This situation makes the use of new tools for nematode identification appealing and necessary.

Molecular approaches to describe, catalog, and identify biological diversity have been increasingly adopted in biodiversity studies, particularly since the inception of the DNA barcoding of life initiative (Hebert et al. 2003a, b; Blaxter 2004). Currently, these methods are being applied to a wide range of taxonomic groups (Floyd et al. 2002; Hebert et al. 2004; Ward et al. 2005; Meyer et al. 2008). The main goal of DNA barcoding is to characterize the biological diversity using a short DNA sequence that will facilitate and expedite taxonomic identification (http://www.barcoding.si.edu/). This initiative has helped in the discovery of new species, many of which have been shown to be morphologically cryptic, thereby considerably improving biodiversity assessment in poorly studied groups such as microscopic meiofauna (Rocha-Olivares et al. 2001; Blaxter et al. 2004; Derycke et al. 2005; Leasi and Todaro 2009). Among marine nematodes, the combination of molecular and morphological approaches has helped to disentangle different species complexes (e.g., Derycke et al. 2008; Fonseca et al. 2008).

DNA barcoding studies of marine nematodes have revealed a good degree of concordance between traditional morphology-based taxonomy and DNA sequences (De Ley et al. 2005; Bhadury et al. 2006, 2008). This is generally reflected in the correspondence between morphologically based Operational Taxonomic Units (OTU) and corresponding Molecular Operational Taxonomic Units (MOTU) (following Blaxter 2004 definitions). The advantages of molecular identification reside in its generally fast and easy implementation, as well as in the potential use of the molecular data for phylogenetic analyses. A partial sequence of the gene coding for the mitochondrial subunit I (COX I) of the cytochrome c oxidase, originally proposed as the standard for DNA barcoding (Hebert et al. 2003b), has so far been unsuccessfully applied in marine nematodes mainly due to the unavailability of PCR primers working across the entire phylum (Blaxter et al. 2005; De Ley et al. 2005; Bhadury et al. 2006; Creer et al. 2010).

Even though no standardized gene for DNA barcoding is available for marine nematodes, an accumulation of nematode 18S or 28S nuclear ribosomal DNA (rDNA) sequences in public databases reflects their usefulness in molecular phylogenetic studies (Nadler 1992; Blaxter et al. 1998; Litvaitis et al. 2000; Nadler et al. 2006; Meldal et al. 2007). These data also reflect that the diverse order Enoplida Filipjev 1929 has been poorly explored compared to other marine nematode taxa by both molecular and classical taxonomy studies. Enoplids are present in nearly all marine sediments and represented by diverse species from several trophic levels (Platt and Warwick 1983). Some of them are predators and play an important role in regulating other nematode populations, mainly in intertidal sandy beaches where they are considered quite common (Greenslade and Nicholas 1991; Nicholas 2002).

The broad application of DNA barcoding to marine nematodes requires first finding a suitable genomic region, or combination of regions, for species identification across a variety of taxa. Second, it requires building a reference database of sequences and morphological vouchers from widespread localities, since most nematode sequences in molecular databases come from NW Europe. Preferably, the DNA barcoding region should have a “barcoding gap,” which refers to a range of sequence divergence between species higher than and non-overlapping with the range of intraspecific divergence (Wiemers and Fiedler 2007). The aim of this study is to assess the marine nematode diversity of the order Enoplida from the coasts of Baja California using an integrative approach. Understanding the levels of concordance between morphological and molecular approaches through integrative taxonomy, as well as addressing the congruence of candidate gene regions for DNA barcoding (i.e., 18S and 28S genes), will help to improve our assessment of marine nematode diversity and their evolutionary relationships.

Materials and methods

Sampling, nematode extraction, and identification

Organisms were sampled in the intertidal zone of sandy beaches with a transparent corer (2 cm Ø × 10 cm) and were immediately fixed in the field with DESS solution (DMSO 20%, 0.25 M disodium EDTA, and saturated with NaCl, pH 8.0; Yoder et al. 2006). Four sites were sampled in the Gulf of California: (1) San Felipe and (2) Santa Clara in the Upper Gulf of California (UGC), (3) Bahía de Los Ángeles close to the “Grandes Islas” region and (4) La Paz Bay at the southwestern end of the Gulf. On the Pacific coast, samples were collected from: (5) San Carlos and (6) Faro beach within Vizcaíno Bay and from (7) Cerritos beach at the southern end of the peninsula (Table 1).

Table 1

Sampling sites

Region	Sitea	Temperature (°C)	Salinity (‰)	Position (lat./long.)	Date	Habitat type
Gulf of California	1. San Felipe (SF)	32	35.5	31°13′41.58″/114°52′48.84″	07/27/07	Pristine, reflexive beach, wide tidal range, gravel sand, small waves
	2. Santa Clara (SC)	30.5	35	31°41′7.98″/114°30′31.74″	07/28/07	Pristine, dissipative beach, wide tidal range, fine sand (308 μm), small waves
	3. Bahía de Los Ángeles (BLA)	NA	NA	28°58′4.20″/113°32′27.00″	08/15/06	Pristine, bay, wide tidal range, fine sand, no wave action
	4. La Paz (LP)	22	36	24° 9′55.62″/110°18′51.00″	11/17/07	Disturbed, bay, fine sand, no wave action
Pacific coast	5. San Carlos (SCA)	NA	NA	29°37′21.72″/115°29′1.62″	05/25/07	Little disturbed, dissipative beach, fine sand, can be affected by big waves
	6. Playa del Faro (PF)	26	36	28°14′21.58″/114° 6′13.62″	05/24/07	Pristine, dissipative beach, median sand, no wave action
	7. Cerritos (CE)	28	35	23°19′44.40″/110°10′32.88″	11/21/06	Little disturbed, dissipative beach, fine-median sand, can be affected by big waves

NA missing data

aAbbreviations in parentheses

Sediment samples were rinsed with tap water on a 63-μm sieve and the meiofaunal community was extracted by floatation with Ludox™ (specific gravity 1.15; Jonge and Bouwman 1977; Somerfield and Warwick 1996). Enoplids were individually picked out with a special brush under a dissection microscope (SZX7 OLYMPUS, 56X) and placed on a slide in a drop of sterile water. These temporary slides were analyzed under a compound microscope (OLYMPUS-BX51) with differential interference contrast. Anatomical details of each specimen were photographed at different magnifications to allow subsequent identification and measurements. Morphological identification of the specimens was based on the latest available keys for Enoplida (Smol and Coomans 2006) and with the help of the database NEMYS (http://nemys.ugent.be/). Morphological vouchers for the specimens (digital microphotographs) are available on-line at the Nematode Tree of Life (NemATOL) web site (http://nematol.unh.edu/). After microscopic observation, which was performed as fast as possible to avoid DNA degradation, nematodes were subject to DNA extraction and PCR.

DNA extraction, amplification, and sequencing

Prior to DNA extraction, each specimen was rinsed (3×) with sterile water to remove traces of DESS. Organisms were then transferred to a sterile slide containing 20 μl of Worm Lysis Buffer (WLB) (50 mM KCl, 10 mM Tris–Cl pH 8.3, 2.5 mM MgCl2, 0.45% NP40, and 0.45% Tween 20 as described in Williams et al. 1992) and 2 μl of proteinase K (10 mg ml−1 stock). Subsequently, organisms were cut into three or more pieces (depending on size) with a sterile scalpel, transferred into 200-μl tube and frozen for 10–30 min at −20°C. Samples were incubated for 1 h at 65°C for protein digestion followed by 10 min at 95°C for proteinase inactivation. Finally, tubes were centrifuged for 1 min at 13,000 rpm. PCR amplifications were performed using 2.5 μl of the extraction supernatant.

Two rDNA genes were partially amplified by PCR: a fragment ca. 350 bp of the 18S gene (small subunit or SSU) and ca. 800 bp spanning the D2–D3 domains of the 28S gene (large subunit or LSU). The SSU gene was amplified using primers MN18F and a degenerate version of 22R (d22R) (GCCTGCTGCCTTCCTTRGA) from Bhadury et al. (2006). The D2D3 region was amplified using primers D2A and D3B from De Ley et al. (2005). Each 25 μl PCR reaction for SSU consisted of 2.5 μl DNA, 10 μl dNTPs (0.5 mM each), 2.5 μl 10× PCR buffer (15 mM MgCl2), 1 μl of each primer (10 μM), 1 U of Taq polymerase (New England Biolabs, 5 U/μl) and ddH2O. The low amplification success of the LSU gene using Taq was partially solved using a high-performance DNA polymerase with the following protocol: 15 μl dNTPs (0.5 mM for each), and 0.75 μl of DyNAzyme™ EXT polymerase (Finnzyme, 1 U/μl).

After verification in a 1.5% agarose gel stained with ethidium bromide (0.5 μg ml−1), PCR products were purified for sequencing using exonuclease and shrimp alkaline phosphatase digestion with EXOSAP-IT (USB Affymetrix, Inc.) following the manufacturer’s protocol. Both rDNA genes were sequenced in both directions with PCR primers using ABI-PRISM® Dye-DeoxyTerminator Big Dye™ v3.1 (Applied Biosystems Inc, CA) with an automatic sequencer Gene Analyzer ABI 3100 (Applied Biosystems Inc, CA).

Data analyses

Morphometric analyses included standard characters used in nematode systematics such as body length (L), body width (W), pharynx length (Ph), nerve ring position (nr), buccal cavity length (bc L), buccal cavity width (bc W), head width (h W), anterior cephalic setae length (acs L), posterior cephalic setae length (pcs L), anal body diameter (abc), tail length (Tail), and shape parameters a (L/W), b (L/Ph), and c (L/Tail). For most taxa with few specimens or with low phenotypic variability, morphometric data (mean and range) were compared with those in the literature for identification. For the three closely related species of Mesacanthion (M. sp1, M. sp2 and M. sp3), statistical analyses were carried out to establish the significance of their morphological distinction, which suggested the existence of several species. For these specimens measurements were used to compute pair-wise Euclidean distances among individuals. Non-metric multidimensional scaling (MDS) was used to assess morphological differentiation on a low dimensional space. Significant differentiation (P < 0.05) among groups was assessed with analysis of similarity (ANOSIM, Clarke and Gorley 2001). Finally, we used analysis of similar percentages (SIMPER) to identify which morphological characters contributed most to the differentiation among groups. Sequential Bonferroni correction was applied to significance levels to adjust for non-independent multiple comparisons (Rice 1989).

DNA sequences were edited with CodonCode Aligner 2.0.1 and subsequently aligned in ClustalX with default parameters (Thompson et al. 1997). Phylogenetic relationships among sequences were estimated with maximum parsimony (MP), maximum likelihood (ML) and neighbor joining (NJ), using heuristic searches in PAUP 4.0b (Swofford 1998). For ML searches, we used the AIC criterion to find the best-fit model of molecular evolution (18S: TVMef+I+G and 28S: GTR+I+G) and its parameters with the programs Modeltest 3.7 (Posada and Buckley 2004; Posada and Crandall 1998) and PAUP 4.0b (Swofford 1998). NJ searches were performed on ML distances computed with the best-fit model of molecular evolution. Non-parametric bootstrap was used to assess branch support (MP and NJ: 1,000 pseudoreplicates, ML: 100 pseudoreplicates).

The incongruence length difference test (ILD test; Farris et al. 1994) was used to determine if tree topologies obtained from both rDNA genes were significantly incongruent. This analysis was carried out using MP heuristic searches and 1,000 bootstrap pseudoreplicates in PAUP 4.0b. Finally, phylogenetic trees were inferred from both genes under MP and ML. For ML, the model of molecular evolution was readjusted for both genes (GTR+I+G) following the AIC criterion in Modeltest 3.7 (Posada and Crandall 1998; Posada and Buckley 2004).

In discussing MOTUs, we adopted a cut-off level of 99.5% sequence similarity (equivalent to no more than 0.5% sequence divergence) among specimens to be assigned to the same MOTU.

Results

Morphological identification

We isolated 139 enoplid specimens from the 7 coastal localities and recognized 22 possible morphological species. From these, 20 were identified to generic level (14 genera, 6 families), one specimen was determined only to family level (Thoracostomopsidae), and another could not be identified beyond the order Enoplida (Table 2). Thoracostomopsidae and Oncholaimidae were the most abundant families with seven and three genera, respectively; Enchelidiidae included two genera, and the rest of the families (Ironidae, Tripyloididae, and Rhabdonemaniidae) only one.

Table 2

Taxonomic list of identified families and genera, number of individuals processed for each morphological species and number of sequences produced for both genes

Order	Family	Genus	Site	Num. (Ind.)	Num. (Seq.)
					18S	28S
Enoplida	Enchilidiidae	Calyptronema	LP	1	0	0
		Pareurystomina	SC	2	0	2
	Ironidae	Trissonchulus	BLA	1	0	0
	Oncholaimidae	Metoncholaimus	PF	8	6	2
		Oncholaimus	LP	9	2	2
		Viscosia	SCA	8	2	2
	Rhabdodemaniidae	Rhabdodemania	SC	9	6	2
	Thoracostomopsidae	Enoploides	SF	20	9	10
		Epacanthion	SC	29	9	3
		Mesacanthion sp1	BLA	10	2	2
		Mesacanthion sp2	SF	15	6	6
		Mesacanthion sp3	CE	11	4	4
		Mesacanthion sp4	PF	1	1	1
		Mesacanthion sp5	SCA	1	1	1
		Mesacanthoides	SC	2	2	1
		Oxyonchus	SC	2	2	1
		Thoracostomopsid	SC	2	2	2
		Trileptium sp1	SC	2	2	1
		Trileptium sp2	SC	3	3	3
		Trileptium sp3	PF	1	1	1
	Tripyloididae	Bathylaimus	SC	1	1	1
	Enoplid	Enoplid	SCA	1	1	1
Total				139	62	48

Site abbreviations are presented in Table 1

Even though organisms identified in this study showed qualitative similarities with species already described in the literature (Trileptium sp2 similar to T. parisetum, Oxyonchus to O. dentatus and O. evelynae, Enoploides to E brunettii and E. longispiculosus, Epacanthion to E. oweni and E. oliffi), morphometric characters suggest that many of our species are new to science (Table 3). In addition, the large unidentified Thoracostomopsidae from Santa Clara (>5 mm in length) differed from all known genera reported in the literature, and was mostly distinguishable by its hemispherical tail, making it a candidate for a novel genus.

Table 3

Morphometric comparison (mean and range in μm) between some of the morphological species identified in this study and species already described in the literature

Species	Sex	L	W	Ph	nr	a	b	c	Tail	V%
Epacanthion sp. from SC (n. sp.), 4 males, 4 females	M	3440 (3330–3710)	70 (70–80)	930 (88–1010)	140	46.0 (44.9–47.5)	3.7 (3.6–3.8)	13.1 (12.2–13.9)	260 (260–270)
	F	3570 (3390–3750)	100 (90–100)	990 (970–1000)	150	37.0 (36.8–37.2)	3.6 (3.5–3.7)	14.4 (13.0–15.8)	250 (240–260)	58.6 (58.0–59.2)
Epacanthion oweni Keppner, 1986, 3 males, 2 females	M	2110 (2050–2170)	51 (50–53)	538 (525–561)	122 (117–128)	41.4 (39.8–43.4)	3.9 (3.9–4.0)	16.5 (16.0–16.9)	128 (127–130)
	F	2580 (2530–2630)	62 (61–63)	593 (589–596)	131 (128–133)	41.6 (41.5–41.7)	4.3 (4.2–4.4)	15.7 (15.7–15.8)	164 (160–167)	59.0
Epacanthion oliffi (Inglis, 1966), 2 males, 3 females	M	1995 (1880–2110)	54.5 (54–55)	615 (590–640)	135.5 (131–140)	36.6 (34.8–38.4)	3.1 (3.2–3.1)	10.2 (10.1–10.3)	195 (182–208)
	F	1943 (1850–2020)	64.3 (58–73)	616.7 (530–680)	155 (138–165)	30 (25.3–32.6)	3.2 (2.9–3.5)	10.9 (9.0–12.1)	182 (138–165)	59.7 (58.4–61.7)
Enoploides sp. from SF, 2 males, 4 females	M	2330 (2170–2490)	70 (60–70)	690 (640–730)	140	35.3 (33.7–36.8)	3.4 (3.4–3.4)	18.6 (17.7–19.6)	130 (120–130)
	F	2470 (2230–2620)	80 (70–80)	700 (620–740)	150 (130–160)	32.5 (31.1–34.9)	3.52 (3.4–3.6)	17.4 (15.7–18.7)	140 (140–150)	65.4 (63.8–66.4)
Enoploides brunettii Gerlach, 1953, 2 males, 1 female	M	2120 (2060–2181)	65.5 (57–74)	460 (450–470)	?	32.5 (29–36)	4.55 (4.5–4.6)	15.35 (14.7–16)	141 (135–147)
	F	2128 (2500–3900)	67 (69–75)	455	?	32	4.7	16	130	59.0 (56.0–59.0)
Enoploides longispiculosus Vitiello, 1967	F	(2400–2900)	(92–112)	?	?	(25–30)	?	?	?	(56–58)
Mesacanthion sp1 from BLA (n. sp.), 4 males, 4 females	M	2000 (1890–2080)	50 (50–60)	420 (380–430)	110 (100–120)	36.6 (33.5–38.6)	4.8 (4.4–5.3)	13.1 (12.4–14.3)	150 (140–170)
	F	1980 (1900–2013)	50 (50–60)	410 (380–460)	110 (100–120)	35.8 (32.3–42.7)	4.8 (4.7–5.0)	12.9 (12.4–13.4)	150 (150–160)	58.3 (57.7–59.0)
Mesacanthion sp2 from SF (n. sp.), 1 male, 1 female, 4 juveniles	M	2870	60	510	130	49.5	5.6	19.0	150
	F	2470	60	460	120	40.6	5.3	18.3	140	52.3
	J	2240 (1970–2540)	50 (50–60)	440 (410–460)	120 (110–130)	44.5 (42.8–46.2)	5.1 (4.7–5.7)	16.5 (14.0–17.9)	140 (120–140)
Mesacanthion sp3 from CE (n. sp.), 4 female	F	1870 (1700–1960)	50	510 (470–550)	110 (110–120)	38.3 (33.9–41.6)	3.6 (3.5–3.9)	14.33 (13.8–15.4)	130 (120–140)	61.4 (60.4–62.0)
Mesacanthion alexandrinus Nicholas, 1993, 4 males, 3 females	M	1830 (1450–2460)	47 (37–59)	509 (430–575)	134 (114–153)	39.2 (33–43)	3.4 (3.1–3.6)	16.5 (14–2)	111.5 (97–110)
	F	2467 (2270–2570)	71 (60–86)	681 (594–730)	164 (151–180)	35.7 (26–42)	3.6 (3.5–3.8)	18 (17–19)	137.3 (136–139)	57.7(55–63)
Oxyonchus sp. from SC, 2 males	M	2560 (2210–2910)	80	750 (720–770)	135 (130–140)	31.9 (28.3–35.5)	3.4 (3.1–3.8)	9.7 (8.8–10.5)	260 (250–280)
Oxyonchus evelynae Nicholas, 2004, 4 males	M	2134 (1964–2397)	76.82 (70.42–81.09)	599 (550–634)	129 (119–146)	46 (38–53)	3.6 (3.3–3.8)	10 (9–12)	225 (195–252)
Oxyonchus dentatus Filipjev, 1927, 1 male	M	4500	100	763	?	45	5.9	15	300
Trileptium sp2 from SC (n. sp.), 1 male, 1 female	M	3310	50	1100	320	66.2	3.0	18.4	180
	F	3500	50	640	160	70	5.5	17.5	200	61.3
Trileptium parisetun Warwick and Platt, 1973, 1 male	M	3800	41	722	?	92.7	5.3	23.1	164.7

Morphometric characters and sites abbreviations were presented in Materials and methods and Table 1, respectively. Data were taken from original descriptions and from Platt and Warwick (1983, in bold). ? denotes missing data

SC Santa Clara, SF San Felipe, BLA Bahía de Los Ángeles, CE Cerritos

Nematodes of the genus Mesacanthion were sorted into five species, based on morphology and molecular data, two in the Gulf of California (M. sp1 and M. sp2) and three in the Pacific coast (M. sp3, M. sp4, and M. sp5). Mesacanthion sp1, M. sp2, and M. sp3 were morphologically similar having a conic-cylindrical tail (n = 8, 6, and 5, respectively). An MDS plot showed a clear segregation of organisms into three morpho-groups suggesting distinct morphological species (Fig. 1). ANOSIM confirmed significant differentiation among them (M. sp1 × M. sp2, P = 0.001 and R = 0.652; M. sp1 × M. sp3, P = 0.001 and R = 0.992; M. sp2 × M. sp3, P = 0.002; R = 0.768; P values significant after sequential Bonferroni correction). Four morphological characters were mainly responsible for this differentiation as revealed by SIMPER. Differences in body ratios a (body length/width), b (body length/pharynx length), and c (body length/tail length) as well as differences in L (body length) accounted for 70–80% of the cumulative differentiation among M. sp1, M. sp2, and M. sp3.

Fig. 1

Non-metric multidimensional scaling (MDS) plot obtained from the morphometric data of Mesacanthion spp. with conic-cylindrical tails. Site abbreviations are in Table 1

Mesacanthion sp1, sp2, and sp3 were morphologically similar to M. alexandrinus in having somatic setae along the body (mainly in the anterior region), a long and slender spicule with a small gubernaculum and one small pocket shape supplement among males. On the other hand, morphometric differences (average and range values) were observed among these species for characters such as body length (M. sp3 < M. sp1 < M. alexandrinus < M. sp2), pharynx length (M. sp1 < M. sp2 < M. sp3 < M. alexandrinus), nerve ring position (M. sp1 < M. sp3 < M. sp2 < M. alexandrinus), and tail length (M. alexandrinus < M. sp3 < M. sp2 < M. sp1, Fig. 2). Based on these comparisons, M. sp2 seems to be more similar to M. alexandrinus than M. sp1; however, they still differ in vulva position (Table 3).

Fig. 2

Whisker plots of mean, standard deviation and range of morphometric variables among three closely related Mesacanthion spp. identified in this study and M. alexandrinus (Nicholas 1993). a Body length, b nerve ring position, c pharynx length, and d tail length

18S and 28S sequences

The DNA of 110 from the 139 analyzed nematodes could be successfully amplified and sequenced (18S = 62 and 28S = 48, Table 2). We identified 19 distinct 18S and 20 28S (D2D3 region) MOTUs. Sequences are available on GenBank under accession numbers: 18S, GU139747–GU139765; 28S, GU139766–GU139785. PCR rate success was different between genes, and was considerably higher (137/139 or 99%) with the short 18S fragment than with the D2D3 region (93/139 or 67%). On the other hand, sequencing success rate was very similar between genes and only a few PCR products could not be successfully sequenced. This was the case of Pareurystomina for the 18S gene and Trissonchulus and Calyptronema for both genes. The alignment of the 18S gene was 324 bp long, including gaps, and the fraction of variable sites was smaller than the fraction of conserved sites (Table 4). On the other hand, the 28S gene alignment was 794 bp long, contained more gaps, and was more polymorphic than the SSU alignment. In this gene, the percentage of variable sites more than doubled the percent of conserved sites (Table 4).

Table 4

Variability and composition of the 18S and 28S rDNA sequences

Gene	Alignment size (bp)	Variable (bp)a	Conserved (bp)	Parsimony informative (bp)b	Nucleotide frequencies (%)
					T(U)	C	A	G
18S	324	138	183	111	27.0	21.5	27.9	23.6
		42.6%	56.5%	34.26%
28S (D2D3 region)	794	545	240	457	24.2	21.4	25.4	29.0
		68.6%	30.2%	57.56%

Alignment including gaps, number of variable, conserved sites and parsimony informative sites (bp and percentage) and nucleotide frequencies (percentage)

aExcluding indels

bIndels were coded as missing data in MP analysis

Comparison of our sequences with those in GenBank revealed high similarities with confamilial and congeneric sequences, and in some cases with sequences from specific taxa. Sometimes 18S and 28S sequences from the same specimen provided matches to different taxa, which likely resulted from the absence of more closely related species in the database (see electronic supplementary material). For instance, Viscosia specimens sampled in San Carlos were 96% similar to the 18S from Viscosia viscosia (DQ394740) with 100% of coverage; 95 and 91% similar to the 28S sequence from Oncholaimus sp. (AF210413) and Viscosia sp. (DQ077751) with 43 and 100% of coverage, respectively.

Morphological and molecular approaches were concordant in species identification. Each morphological species with replicate specimens possessed a unique rDNA sequence in both genes resulting in a single MOTU. Although not all the specimens were sequenced, most of the species had two sequences for intraspecific comparison, at least for one of the genes (Table 2). MOTU divergences increased with decreasing relatedness in both genes; however, divergence was consistently higher for the 28S rDNA (Fig. 3).

Fig. 3

Frequency distribution of percent sequence divergence in Enoplida for 18S (white bars) and 28S (black bars) rDNA. Number of pair-wise comparisons: 18S: 1770, 28S: 1035

We found morphological variation among specimens of Enoploides and Mesacanthion sp3 potentially interpretable as the presence of more than one species. In Enoploides, a male showed features (head shape, jaws, and body size) very different from the rest of congeneric males and females. Mesacanthion sp3 specimens, on the other hand, were entirely females and juveniles with very different morphological features (mainly body shape and head, electronic supplementary material); nevertheless, they shared identical DNA sequences of both genes. The same was true for the specimens of Enoploides. Consequently, this morphological variation was not interpreted as indicative of additional taxa.

Phylogenetic analyses

The three methods of phylogenetic reconstruction (ML, MP, and NJ) produced highly congruent 18S gene trees, featuring three major clades (Fig. 4a). A highly supported clade (bootstrap > 90%) consisted of all the Thoracostomopsidae sequences; a second grouped Bathylaimus and (non-identified) Enoplida n.i. (bootstrap 100%); a third moderately supported clade grouped the three oncholaimid genera (Viscosia, Metoncholaimus, and Oncholaimus) with Rhabdodemania (Rhabdonemaniidae) (Fig. 4a). Thoracostomopsidae and Oncholaimidae were monophyletic in all trees, and the latter received the highest bootstrap support in all reconstructions. Whereas the three oncholaimid genera were completely resolved, some relationships among the more numerous Thoracostomopsid genera were not. The five sequences of Mesacanthion spp. were not monophyletic. Whereas the closely related species M. sp1 and M. sp2 were grouped with a maximum support, Epacanthion and M. sp3 were joined to the same clade with modest support values (Fig. 4a).

Fig. 4

Maximum likelihood phylogenetic reconstructions of Enoplida sequences from Baja California. a 18S gene, b 28 S gene. Bootstrap values (only >50%) for the different methods (MP, ML and NJ, respectively) are shown on branches (– bootstrap <50% or no value at all, * branch absent). Families are indicated by vertical bars

On the other hand, 28S gene trees were less consistent. The main difference among the inference methods was the position of the enchilidiid genus Pareurystomina (Fig. 4b, in bold). In the MP reconstruction, this genus was sister to Oncholaimidae and relationships were completely resolved with families reciprocally monophyletic. The relationship of Viscosia as sister to Oncholaimus and Metoncholaimus was moderately supported (60% bootstrap). In the ML and NJ trees, Oncholaimidae and Enchilidiidae were paraphyletic, although relationships within these families were poorly resolved. However, the taxa comprising both families grouped in a monophyletic Oncholaimoidea (Fig. 4b).

Since only one Enchilidiid was present in our samples, an additional sequence was obtained from GenBank (28S: Calyptronema maxweberi, AF210399) to address relationships between Enchilidiidae and Oncholaimidae. The reciprocal monophyly of the families was not supported by our data: MP produced a polytomy joining both enchilidiids (Pareurystomina and C. maxweberi) with the monophyletic oncholaimid sequences; ML produced a paraphyletic Enchilidiid pair (Fig. 5b). The monophyly of the superfamily Oncholaimoidea (Enchilidiidae + Oncholaimidae) in these analyses was consistently and strongly supported (Fig. 5).

Fig. 5

a Maximum parsimony and b maximum likelihood phylogenetic reconstructions of Enoplida 28S rDNA sequences from Baja California including Calyptronema maxweberi sequence from GenBank (AF210399.1). Bootstrap values (only >50%) are shown on branches

The phylogenetic signal contained in each rDNA gene produced topological differences in the phylogenies. For instance, the genus Rhabdodemania was sister to Enoplida n.i. + Bathylaimus in the 28S tree but sister to all Oncholaimidae in the 18S reconstruction; both nodes were well supported by at least two of the three methods (MP, ML, NJ). Other differences were found among Thoracostomopsid genera, particularly in the deeper nodes leading to their common ancestor. Based on the number of nodes with high support values, both genes appear to have comparable resolving power but for different taxa.

In spite of these topological differences, the 18S and 28S trees were congruent (ILD test, P = 0.28), suggesting that differences were the result of poorly resolved relationships in both data sets. Therefore, we proceeded to merge the sequence data for a combined analysis (1,123 characters, including gaps). MP and ML (evolution model readjusted to GTR+G+I) produced very similar trees with three very well-supported main clades: (1) Thoracostomopsidae, (2) Oncholaimidae and (3) Rhabdodemania, Bathylaimus and Enoplida n.i. (Fig. 6); of these only the latter received less than 90% bootstrap support in the ML reconstruction (Fig. 6b). Within Thoracostomopsidae, ML produced a larger number of better supported clades, including the monophyly of the three closely related Mesacanthion spp., which were paraphyletic in the MP consensus tree (Fig. 6a). Trileptium sequences were consistently paraphyletic in all phylogenetic reconstructions. Only the monophyly of sympatric Trileptium sp1 and sp2 was recovered with high support values in all analyses, whereas Trileptium sp3 from Playa del Faro was paraphyletic with Mesacanthion sp4 and M. sp5 in both trees (Fig. 6).

Fig. 6

Phylogenetic reconstruction using 1,123 bp from 18S and 28S rDNA genes combined of Enoplida nematodes from Baja California. a MP consensus tree of five most parsimonious reconstructions. b ML tree based on GTR+G+I model of sequence evolution. Bootstrap values (only >50%) are shown on branches

Discussion

The marine nematode communities of the Gulf of California are poorly known and data are available from only two localities in the UGC (Mundo-Ocampo et al. 2007; Holovachov et al. 2008), and there is no information available from the Pacific coast of Baja California. Thus, this study focused on enoplids reveals novel insights about the distribution of marine nematodes in other areas of the Gulf of California as well as the first data from the Pacific coast of Baja California, Mexico.

Integration of morphological and molecular approaches for identification

Morphological identification of the enoplids in this study was supported by the molecular data, as shown by the congruence and high similarity between our sequences and those available in GenBank. Discrepancies between our morphological identifications and the closest matches in the molecular database are most likely the result of the limited taxonomic coverage of rDNA sequences available in GenBank relative to the vast diversity of marine nematodes (Lambshead and Boucher 2003). The limited sequence availability and the unavailability of morphological cross-referenced vouchers in GenBank preclude using these Blast analyses as a bona fide molecular identification tool. Pending the sampling of additional molecular data of closely related taxa, they are nevertheless consistent with our morphological determinations.

In addition, molecular data also revealed that some conspicuous morphological variation between congeneric specimens could be ontogenetic (M. sp3, differences between juvenile and adult) or due to phenotypic plasticity (Enoploides, differences between adults). Ontogenetic variation is an important issue in nematode identification, since most diagnostic characters relate to adult (often male) genitalia absent in subadults, which may be the only life-stage sampled (De Ley 2000; Bhadury et al. 2006). On the other hand, marine nematodes seem to be phenotypically plastic, which is a source of taxonomic uncertainty and may mislead and bias diversity studies based exclusively on morphology (Nadler 2002; Derycke et al. 2008). Additional molecular data from more polymorphic genes would be required to test the hypothesis that rDNA sequences were insufficiently variable to reflect interspecific divergence if the observed phenotypic variation relates to differences between species instead of intraspecific plasticity.

Surprisingly, two of the Mesacanthion spp. (M. sp1 and M. sp2) found in the Gulf of California presented morphological features very similar to M. alexandrinus described by Nicholas (1993) from a freshwater environment in Australia. Despite their qualitative resemblance, the contrasting habitats where they were found (freshwater vs. marine) and the levels of phenotypic differentiation suggest that Mesacanthion spp. from this study may represent species new to science. Based on the integration of morphological and molecular data, we also showed that the three closely related Mesacanthion spp. (M. sp1, M. sp2 and M. sp3) are different from each other and therefore should be treated as different species. The unidentified Thoracostomopsidae differed from all genera known and reported in the literature. We could find no resemblance to any genus described for the entire order Enoplida, even to those described from freshwater environments (Smol and Coomans 2006), suggesting that this nematode may represent a new genus of Thoracostomopsidae. A detailed morphological description of these specimens is beyond the scope of this paper and will be presented elsewhere.

Morphological and molecular approaches were congruent in addressing nematode species identification since (1) each species identified based on morphology presented a unique DNA sequence (MOTU), and (2) sequences were phylogenetically concordant with taxonomy, for most part. Moreover, the combination of both approaches showed that natural variability (ontogenetic development, sibling and cryptic species) could bias biodiversity evaluation, over- or underestimating the number of species. An integrative taxonomic approach is the best strategy for marine nematode identification as previously suggested (Coomans 2002; De Ley 2006; Bhadury et al. 2008; Derycke et al. 2008).

Relative merit of 18S and 28S genes for DNA barcoding

Both genes, 18S and 28S, produced an equal number of unique sequences (MOTU diversity) showing the same capability of addressing nematode species identification. The smaller 18S fragment showed a considerably higher amplification success; however, the lower polymorphism may limit its phylogenetic usefulness. On the other hand, divergence was much higher in the D2D3 domain of the 28S gene, regardless of taxonomic level (i.e. families or genera), making it suitable for both DNA barcoding and phylogenetic reconstruction of marine nematodes (De Ley et al. 2005; Derycke et al. 2008). This increased divergence would be particularly valuable in resolving closely related species, as in the case of Mesacanthion (M. sp1 and M. sp2, 18S = 1.06% and 28S = 12.09%) and Trileptium (T. sp1 and T. sp2, 18S = 3.87% and 28S = 12.97%), and to detect possible cryptic species, where morphological differences are sometimes impossible to diagnose (De Ley et al. 1999; Fonseca et al. 2008). However, 28S had limited success in PCR amplifications (67%). A gene used for DNA barcoding should be easily reproducible across the entire Nematoda phylum. The value of the D2D3 region as a potential barcoding gene has been shown in studies of plant parasitic and free-living nematodes (De Ley et al. 1999; Subbotin et al. 2005; Subbotin et al. 2007). De Ley et al. (2005) have shown a high success rate in PCR amplification between different groups of nematodes; consequently, increased success for marine nematodes may require additional PCR optimization or the design of new nested primers encompassing diagnostic regions (Fonseca et al. 2008).

Phylogeny

In general, phylogenetic inferences from both genes were not intrinsically different among methods (MP, ML and NJ). In fact, tree topologies for the 18S gene were completely congruent. Differences among 28S reconstructions involved only the position of the Enchilidiid Pareurystomina (Fig. 4b). Inclusion of an additional Enchilidiid sequence (28S: Calyptronema maxweberi, AF210399.1) did not help resolving the reciprocal monophyly of Enchilidiidae (Pareurystomina and C. maxweberi) and Oncholaimidae (Viscosia, Oncholaimus, and Metoncholaimus) (Fig. 4). Additional taxa from both families may be needed to resolve this node, should the lack of resolution be a result of insufficient taxon sampling.

This study also revealed contrasting phylogenetic signals in 18S and 28S genes. Although the ILD test did not detect significant differences between tree topologies (P = 0.28), the genus Rhabdodemania showed a controversial position in the 18S reconstruction. In morphology-based and other molecular phylogenies (e.g., Litvaitis et al. 2000), Rhabdodemania (Rhabdodemaniidae) and Bathylaimus (Tripyloididae) are grouped in the suborder Tripyloidina, in concert with the topology obtained with the D2D3 fragment (Fig. 4).

Our results underline the need for a combination of morphological and molecular approaches to expedite our understanding about marine nematode taxonomy, biogeography, dispersal capacity and gene flow among populations. These approaches will prove invaluable to have a fresh understanding of real levels of cosmopolitanism in marine nematodes, as it has been extensively reported in the literature (Bhadury et al. 2008; Derycke et al. 2005; Fonseca et al. 2006; Heip et al. 1985). Finally, an integrative approach will aid in the detection of cryptic species, which are common among meiofaunal groups, thereby improving the assessment of marine nematode diversity and contributing to a more robust nematode taxonomy (Derycke et al. 2007; Rocha-Olivares et al. 2001; Todaro et al. 1996).

Supplementary material (DOC 184 kb)