Mitogenomics reveals high synteny and long evolutionary histories of sympatric cryptic nematode species

Abstract Species with seemingly identical morphology but with distinct genetic differences are abundant in the marine environment and frequently co‐occur in the same habitat. Such cryptic species are typically delineated using a limited number of mitochondrial and/or nuclear marker genes, which do not yield information on gene order and gene content of the genomes under consideration. We used next‐generation sequencing to study the composition of the mitochondrial genomes of four sympatrically distributed cryptic species of the Litoditis marina species complex (PmI, PmII, PmIII, and PmIV). The ecology, biology, and natural occurrence of these four species are well known, but the evolutionary processes behind this cryptic speciation remain largely unknown. The gene order of the mitochondrial genomes of the four species was conserved, but differences in genome length, gene length, and codon usage were observed. The atp8 gene was lacking in all four species. Phylogenetic analyses confirm that PmI and PmIV are sister species and that PmIII diverged earliest. The most recent common ancestor of the four cryptic species was estimated to have diverged 16 MYA. Synonymous mutations outnumbered nonsynonymous changes in all protein‐encoding genes, with the Complex IV genes (coxI‐III) experiencing the strongest purifying selection. Our mitogenomic results show that morphologically similar species can have long evolutionary histories and that PmIII has several differences in genetic makeup compared to the three other species, which may explain why it is better adapted to higher temperatures than the other species.


Introduction
Species with similar morphologies but distinct genetic differences (= cryptic species) abound in all animal taxa and in all geographical areas and indicate that morphological stasis is a constant evolutionary phenomenon (Pfenninger and Schwenk 2007). Allopatric speciation has been the basic mechanism to explain species formation (Coyne and Orr 2004), even in marine environments where barriers to gene flow are less obvious, population sizes are large, and dispersal is substantial (Knowlton 1993;Palumbi 1994). Cryptic species have long been recognized in marine environments and are often found together (Knowlton 1993;Stuart et al. 2006). This may point to secondary contact of species that speciated allopatrically, or to the occurrence of species formation in sympatry (Via 2001). In the latter scenario, monophyletic sister clades are expected to occur in sympatry (Coyne and Orr 2004). The discovery of cryptic species is generally a side product of population genetic and phylogeographic studies, resulting in a delineation that is often based on a phylogeny generated with a limited number of loci. As such, alterations in gene order, content, and size that may occur during or after the speciation event remain undetected. Mitochondrial (mt) genomes of animals are particularly interesting to investigate evolutionary relationships between closely related species because they evolve faster than nuclear genomes, their small size, simple structure, absence of recombination, and high variability (Castellana et al. 2011).
The metazoan mt genome is typically characterized by a circular double-stranded moleculebut exceptions have been observed, for example, in the nematode Globodera pallida where six small circular molecules were found (Armstrong et al. 2000;Gibson et al. 2007). It encodes for 2 rRNAs, 22 tRNAs, and 13 proteins that are essential for mitochondrial functioning (Boore 1999). The proteins are typically under strong purifying selection (Meiklejohn et al. 2007;Castellana et al. 2011) which prevents the accumulation of mutations that may alter the functional protein products of important subunits of the respiratory chain. Variability in metazoan mitochondrial gene content rarely involves protein-coding and rRNA genes, and is mostly attributed to tRNA genes (Gissi et al. 2008). Phylogenetic relationships based on mitochondrial DNA can be blurred by the presence of Wolbachia which can cause introgression and/or selective sweeps (Hurst and Jiggins 2005). Nevertheless, a mitogenomic approach has yielded novel insights in phylogenetic relationships of many organisms, including nematodes (Sun et al. 2014;Zasada et al. 2014).
Nematodes form one of the most successful animal phyla in terms of species diversity and habitat exploitation (Bongers and Ferris 1999). Substantial cryptic diversity has been observed in parasitic (de Le on and Nadler 2010), freshwater (Ristau et al. 2013), and marine species (Derycke et al. 2013). The genetic differences underlying this cryptic speciation have however been restricted to the use of a few loci (the mitochondrial cytochrome oxidase c subunit 1 (coxI) gene, and fragments of the nuclear ribosomal units). Population genetic and phylogeographic studies have demonstrated the presence of at least 10 cryptic species within the Litoditis marina morphospecies formerly known as Rhabditis (Pellioditis) marina (Bastian 1865)in the northeast Atlantic (Derycke et al. 2008a,b). Increasing sympatry was observed when species were more distantly related, suggesting a predominant mode of allopatric speciation (Derycke et al. 2008a,b). Four of these species (PmI, PmII, PmIII, and PmIV) cooccur on decomposing macroalgae in the littoral zone of coastal and estuarine environments in Belgium and the Netherlands (Derycke et al. 2005). The four species lack diagnostic morphological features and are reproductively isolated . Temporal fluctuations in their distribution have been observed (Derycke et al. 2006). In addition, differences in competitive ability (De Meester et al. 2011), timing of dispersal (De Meester et al. 2012), and microbiomes (Derycke et al. submitted) have been demonstrated. Although their life histories and ecological interactions are beginning to be unraveled, the evolutionary history of this species complex remains lar-gely unknown. Phylogenetic relationships based on the nuclear ribosomal ITS and 28S regions showed that PmI and PmIV are sister taxa, which are more closely related to PmII and more distantly related to PmIII (Derycke et al. 2005(Derycke et al. , 2008a. Although the four species were recovered in the mitochondrial COI gene phylogeny, deeper relationships were not supported (Derycke et al. 2008a,b). In this study, we investigated the whole mitochondrial genomes of these four sympatrically distributed cryptic species and specifically aimed at identifying differences in genomic architecture (size, gene order, codon usage) between the four species. In general, high synteny and chromosomal organization are observed between rhabditid species (Nigon and Dougherty 1949;Hillier et al. 2007;Bik et al. 2012;Sun et al. 2014), although Caenorhabditis briggsae, the cryptic sister species of C. elegans, shows some unusual nad5 pseudogenes and associated heteroplasmic deletion events that suggest a dynamic evolution between C. elegans and C. briggsae (Howe and Denver 2008;Phillips et al. 2015). Second, we investigated phylogenetic relationships using all mitochondrial protein-coding genes. Based on the phylogenetic relationships established from two nuclear loci, we expected PmIII to be most distantly related to the other species and to recover PmI and PmIV as sister taxa. Third, we estimated the timing of divergence between the four species. Phylogeographic results have suggested that this timing was before the last glacial maxima in the Pleistocene (Derycke et al. 2008a,b), but finer time estimates are lacking. Finally, we explored selective pressures on the proteincoding genes (PCGs), and we expected to observe purifying selection because of the important functions of mitochondrial protein-coding genes in the respiratory chain.

Nematode cultures
Four monospecific cultures of Litoditis marina were kept in the lab under stable conditions at 18°C. Species PmI, PmII, and PmIII were isolated from Fucus fragments from Paulina (51°21 0 N, 3°49 0 E), a saltmarsh in the polyhaline area of the Westerschelde estuary (the Netherlands). Species PmIV was isolated from Lake Grevelingen (51°44 0 N, 3°57 0 E), a marine lake in the Netherlands. Each culture was started from a single female and was maintained in the lab for many generations before the experiment. Nematodes were cultured on marine agar plates (25 psu, 4/1 bacto/nutrient agar in a final concentration of 1%) to which 2 mL of Tris-HCl (pH 8) was added to buffer the pH of the plates, and which were seeded with Escherichia coli K12 as a food source. All females were morphologically identified using a stereomicroscope and the species description of Rhabditis (Pellioditis) marina (Bastian, 1865) as outlined in (Inglis and Coles 1961). Species identity of the cultures was checked via qPCR . After successful establishment of the cultures, DNA was extracted from a single nematode by transferring the worm to a 0.5-ml Eppendorf tube containing 20 lL worm lysis buffer (50 mmolÁL KCl, 10 mmolÁL Tris pH 8.3, 2.5 mmolÁL MgCl2, 0.45% NP 40, 0.45% Tween-20). Tubes were frozen at À20°C to disrupt cuticula and cell membranes, after which 1 lL of proteinase K (10 mgÁml À1 ) was added. Lysis was performed by incubating the tubes for 1 h at 65°C and 10 min at 95°C. DNA samples were centrifuged for 1 min at 20,800 g and were then used as template for qPCR.
Assignment of nematodes to one of the four cryptic L. marina species were performed with qPCR using species-specific primers located in the ribosomal internal transcribed spacer (ITS) region . The qPCRs were prepared in 10 lL volume containing 5 lL SensiMix SYBR No-ROX One-Step (29) solution, 3 lL of each primer (final concentrations of 1 lmolÁL for Pm I and Pm III, 500 nmolÁL for Pm II, and 200 nmolÁL for Pm IV ), 1 lL PCR-grade water, and 1 lL of DNA template. The thermal cycling protocol consisted of an initial denaturation for 10 min at 95°C followed by 40 cycles of denaturation for 10 s at 95°C, annealing for 20 s at 60°C, and extension for 20 s at 72°C. All specimens were analyzed with all four primers, and the primer set yielding a positive signal was used to identify the specimen as PmI, PmII, PmIII, or PmIV.

DNA extraction and next-generation sequencing
Nematodes were removed from the agar dishes by sucrose washing (using sucrose in a final concentration of 40%) and washed four times in artificial seawater (ASW). DNA was extracted from several hundred nematodes of each species using the CTAB protocol as described in . DNA concentrations were measured with a Nanodrop ND2000, and quality of the DNA was checked by gel electrophoresis before sending the samples to the Hubbard Center for Genome Studies (University of New Hampshire, USA). Sequencing libraries for the four samples were generated following the "Low Throughput Sample Protocol" for Illumina TruSeq DNA libraries. Approximately 1 lg of DNA was sheared by ultrasonification using a Covaris M22 to a target size of 500 bp. The overhangs resulting from the shearing were converted to blunt ends using the End Repair Mix of the Illumina kit. The resulting samples were cleaned with AMPure Beads, after which a single "A" nucleotide was added to the 3 0 ends of the fragments using the A-Tailing mix of the Illumina kit. For each individual species, a unique adaptor/index was ligated to the DNA fragments after which the samples were cleaned with AMPure Beads. The ligation products were then loaded on a 2% agarose gel to remove unligated adapters and adapter dimers and to size-select the fragments at 500 base pairs. Samples were then purified with the MinElute Gel Extraction Kit (Qiagen, Benelux, Antwerp, Belgium). DNA fragments with adapters on both ends were then enriched by PCR with primers that anneal to the ends of the adapters and using an initial denaturation of 98°C for 30 s, 10 cycles of 98°C for 10 s, 60°C for 30 s and 72°C for 30 s, and a final extension of 5 min at 72°C. PCR products were cleaned with AMPure Beads and analyzed for size and concentration on an Agilent Bioanalyzer. The four libraries were quantified through qPCR (Kapa Biosystems), pooled, and loaded on two lanes of a single rapid-run flow cell for paired-end sequencing (2*150 bp) on the Illumina HiSeq 2500.

Data analysis
Assembly and annotation of the mtDNA A de novo assembly was generated with the CLC main workbench software (www.clcbio.com). As we were mainly interested in the mitochondrial DNA, we filtered the results table containing the contigs from the assembly, with filters set on consensus length >230 and total read count >5. Litoditis marina is a rhabditid nematode and belongs to the same family as the model organism Caenorhabditis elegans for which a completely annotated genome is available. The mitochondrial genome of C. elegans was downloaded from Genbank (accession number NC_001328.1). A local BLASTn search was performed with the mtDNA of C. elegans as query, using default settings except for the e-value, which was set at 1e-10, and the number of threads, which was set at 50. From the retrieved contigs, those in the size range of mtDNA (13,000 bp) were retained for further analyses. GC content for the mtDNA was calculated in R 3.0.1 (http:// www.r-project.org/) using the seqinR package (Charif and Lobry 2007).
Open reading frames (ORFs) in the contigs containing the mitochondrial genomes were searched using CLC. Annotation of the ORFs was performed using a BLASTn search against Genbank and with e-value set at 1e-10 and specified for Metazoa. CLC detected 11 ORFs and the lacking protein-coding gene, nad4L, was found by a specific query using the protein-coding gene sequence of C. elegans (extracted from the complete mt genome of C. elegans accession number NC_001328.1) against our mtDNA contig. Once the 12 expected ORFs were anno- tated, they were verified by checking the start and stop codon positions and by translating the sequence to protein. For all genes, pairwise distances were calculated (p-distribution model, pairwise deletion) between all four species in MEGA6 (Tamura et al. 2013). Codon usage was calculated for the four species in MEGA6 and the Relative Synonymous Codon Usage (RSCU) (Sharp et al. 1988) values were used to compare between species. Gene density of the mitochondrial genomes was calculated by dividing the number of genes by the genome size. The genome size was determined with exclusion of the ATrich region because assembly algorithms generally cannot assemble completely through the AT-rich regions and there may be numerous small contigs for this region that are not included in the main contig. Consequently, variability in length may be caused by an inadequate assembly rather than by actual differences in genome size between the species.
tRNA detection was performed in several steps. First, all the tRNA sequences were extracted from the C. elegans complete mt genome (accession number NC_001328.1) and then queried via BLASTn against the mtDNA of the four cryptic species to find their positions. The secondary structure of these tRNAs was determined in CLC, and we checked that the anticodon was located in the loop (further referred to as manual procedure). Second, the mtDNA was uploaded in tRNAscan-SE 1.21 (Lowe and Eddy 1997) to verify the detected tRNA, with source settings on Nematode Mito. Only the two genes for tRNA Ser were not confirmed with tRNAscan-SE. Finally, the mtDNA sequence was also imported in ARWEN (Laslett and Canb€ ack 2008), another tRNA detection tool, in a further attempt to identify the tRNA Ser not identified with the previous algorithms. With ARWEN one of the two tRNA Ser could be identified.
The annotated mitochondrial genomes have been submitted to Genbank under Accession numbers KR815450 for PmI, KR815451 for PmII, KR815452 for PmIII, and KR815453 for PmIV.

Mode of selection on PCG
Synonymous (dS) and nonsynonymous (dN) substitution rates were calculated in MEGA6 (Tamura et al. 2013), using the Nei-Gojobori method and 500 bootstraps. When dN < dS, purifying selection eliminates new variants. When dN > dS, new variants are selected, and positive selection is acting. Under neutral evolution, dS is equal to dN. Afterward, a Z-test of selection was performed as implemented in MEGA6, with the hypothesis tested being purifying selection, as dN < dS, and following settings were used: 500 bootstraps and the Nei-Gojobori method for synonymous-nonsynonymous substitution type.

Phylogenetic relationships and timing of divergence
All PCG sequences of the four species were extracted, translated to amino acids and individually aligned in CLC, including C. elegans and C. briggsae homologs (downloaded from GenBank, NC_001328.1 and NC_009885.1). The aligned amino acid sequences of the PCGs were then untranslated to nucleotide sequences and concatenated into one fasta file, which was used to determine phylogenetic relationships in MEGA6 (Tamura et al. 2013). Neighbor-joining settings included 500 bootstrap replications, the uncorrected p-distance model and pairwise deletion of gaps, Maximum likelihood settings were GTR + G + I model and 500 bootstrap replications, and Maximum Parsimony settings were 500 bootstrap replications. Default settings were used for the other parameters. The evolutionary model that best fitted our data was obtained using default settings and treating gaps/missing data with partial deletion in MEGA6.
Additionally, the most recent common ancestor (MRCA) was calculated with BEAST v2.1.1, containing the BEAST, BEAUti, and TreeAnnotator programs (Bouckaert et al. 2014). Only the PCGs were used to make the alignment, which was then imported in BEAUti as a nexus file to be able to set the model parameters for BEAST. At the site models tab, the Gamma Category Count was set at 4, the Shape and Substitution rate was estimated and the GTR model was chosen as the substitution model (see above) with an empirical frequency. A normal relaxed clock model was selected, and the priors were set at the calibrated Yule model, and an extra criterion was added to define the calibration node. Divergence time of C. elegans and C. briggsae has been estimated at 18 MYA (Cutter 2008), and both species were set to be monophyletic within Caenorhabditis from this time point by defining a normal distribution with parameters (18, 0.5). The Markov Chain Monte Carlo settings were set at a chain length of 1,000,000, a trace log of 200 and screen log of 1000. TreeAnnotator was used to obtain an estimate of the phylogenetic tree. A 1% burnin was specified (being 50), the posterior probability limit was set at zero, and mean heights were chosen for the nodes. The generated tree file was visualized in FigTree v1.4.0. A second calibration point with the divergence between Chromadorea/Enoplea between 532 and 383 MYA (Rota-Stabelli et al. 2013) was used for comparison.

Wolbachia detection
The assembled contigs of the four species were searched for the presence of the endosymbiotic bacteria Wolbachia spp. Five characteristic genes (coxA, gatB, hcpA, ftsZ, and fbpA), the wsp gene, and the 16s rDNA gene (Doudoumis et al. 2012) of Wolbachia endosymbionts that have been found in the parasitic nematodes Brugia malayi and Onchocerca spp. were downloaded from Genbank (respective accession numbers for coxA are DQ842273 and FJ390245; for gatB DQ842421 and JX075229; for hcpA DQ842384; for ftsZ AY583309 and AJ276501; for fbpA DQ842347 and JX075225; for Wsp AY527201 and AY095210 and CU062443; and for the 16S rDNA AF051145 and AF172401). These sequences were queried against our contig databases containing all assembled contigs using BLASTn. A match was considered when the % identity score was higher than 80.

Comparison of the Litoditis mtDNA genomes
The mtDNA of each species was located on a single contig and was 13,766, 13,855, 14,481, and 13,909 bp long, for PmI, PmII, PmIII, and PmIV, respectively ( Table 2). The mitochondrial genomes encode for 2 rRNAs, 22 tRNAs, and 12 proteins, which were positioned in the same order in each of the four L. marina species (Fig. 1). The assembled AT-rich region showed considerable length variability between the four species and was considerably longer in PmIII ( Table 2). The mitochondrial genome size excluding the AT-rich region showed small differences in length between the four species, with the mt genome of PmIII still being 140-162 bp longer than that of the other species (Table 2). GC content was very similar between the four species and varied between 20.4 and 21.3%. PmIII had a somewhat lower percentage of coding information compared to the three other species (97.5% vs >98.2%) which was caused by the additional DNA sequence found in the intergenic region located between tRNA Met and tRNA Asp .

Protein-coding genes (PCGs)
The gene order was the same in the four species (Fig. 1) and introns were absent in the PCG. Gene sizes of PmI and PmIV were identical, while Atp6 and nad4L were one amino acid (AA) longer and shorter, respectively in PmII. In PmIII, nad5 and coxIII were 1 AA longer than in the other three species and nad3 was 13 AA longer ( Table 3). The ATP synthase F0 subunit 8 (atp8) gene was lacking in the mitochondrial genomes of all four cryptic species.
Codon usage was very similar for the four species and the same preference toward one codon over other codons encoding for the same AA was present (Table 4). Only three codons were never used: CUC(L), CGC(R), and CGG(R). The ATT, ATA, and TTG start codons were used by all four species, while the ATG start codon was observed only once for the nad5 gene in PmIII. The TAA stop codon was used for termination of all PCGs in PmII and PmIII, while the TAG stop codon was used once by PmI and PmIV for termination of the coxIII gene. Differences in the use of start and stop codons were observed between the four species in three genes: nad5, coxII, and coxIII. For nad5, PmIII used the start codon ATG, while the three other species used the ATT codon; for coxII, PmI and PmII used the start codon ATT, while PmIII and PmIV used the start codon ATA. For coxIII, PmI and PmIV had the same start and stop codon (ATT/ TAG), PmII and PmIII used the same stop codon (TAA), but had different start codons, respectively ATT and TTG. For coxII, PmIII and PmIV used a different start codon (ATA), while PmI and PmII used ATT as start codon. P-distance values based on all PCGs ranged between 6.1 and 10.5% within the L. marina species complex ( Table 5). As a comparison, the P-distance between C. elegans and C. briggsae was 13.9%, and values between Litoditis and Caenorhabditis ranged between 17.1 and 17.7%. The highest variability was found in the nad genes (Table 5).

Ribosomal and transfer RNA
The rRNAs were located on the same position for all four species: rrnS between tRNA Glu and tRNA Ser(Undef.) and rrnL between tRNA His and nad3. Small differences in length were observed between the four species, ranging between 639 and 697 for rrnS and between 956 and 960 bp for rrnL (Table 3). The four mt genomes contained 22 tRNA genes. None of them had the conventional cloverleaf structure and showed a reduced TwC stem-loop region. In all four species tRNA Ser(TCT) lacked a D-stem and had a reduced TwC stem, while only tRNA Ser(TGA) of PmIII had a D-stem and a TwC stem. For the other three species, the tRNA Ser(TGA) could not be confirmed and its possible location was determined based on a BLASTn search with C. elegans as a query. This tRNA is named tRNA Ser(Undef.) , due to the unknown anticodon. Some tRNAs had overlap of 1-30 nt with adjacent genes. In all four species the largest overlap was observed for tRNA Thr with 30 nucleotides, followed by an overlap of 3 nt for tRNA Tyr . PmI, PmII, and PmIV also had an overlap of tRNA Ile and tRNA Trp in common with 1 nt. Overlap in tRNA Ser(Undef.) with 1 nt was found in PmI and PmIV. PmIII had an overlap of 3 nt and 10 nt for tRNA Leu(TAG) and tRNA Ser(TGA) , respectively.

Noncoding regions
The AT-rich region is a highly variable region with 310 variable positions in a 958-bp-long alignment. It is located between tRNA Pro and tRNA Ala in all four species, but showed pronounced differences in length for species PmIII compared to the three other species (PmI: 378 bp; PmII: 479 bp; PmIII: 953 bp; and PmIV: 543 bp). An intergenic region located between coxI and nad4 in all four species is also highly variable (45 positions/ 149 bp). The length of the intergenic region was similar for PmI, PmII, and PmIV, being respectively 147, 149, and 146 bp long. PmIII was different in that it had two intergenic regions, one located between coxI and nad4 like the other three species but considerably shorter (85 bp long), and one (206 bp long) located between tRNA Met and tRNA Asp .

Phylogenetic relationships and timing of divergence
The PCG alignment was 10380 bp long. The Litoditis marina complex contained 1631 variable and 261 parsimony-informative positions. Maximum likelihood, neighbor-joining, and maximum parsimony trees gave the same topology, and each branch was supported by bootstrap values of 100. PmI and PmIV were sister taxa, while  PmIII was most distantly related to the three other species (Fig. 2). The time tree calculated in BEAST (Fig. 3) showed that the most recent common ancestor for the Litoditis marina species complex was situated at 16 million years ago (MYA). PmII diverged 10.5 MYA from PmI and PmIV, while the latter two diverged 6.5 MYA. Time calculations were the same when using one or two calibration points.

Wolbachia detection
A BLASTn search was performed with the Wolbachia sequences as queries against the complete genomes of the four cryptic species, containing both nuclear and mtDNA. Of the five conserved genes of the Multilocus Sequence Typing system (MLST), two genes showed no significant match, while the three other genes were found with a low % identity (<75%). Also the Wolbachia surface protein (wsp) and the 16S RNA gene were not detected in the sequence data of the four species. These results suggest that Wolbachia was not present in the genome assemblies of these species.

Mode of selection on PCGs
The ratio of nonsynonymous (dN) to synonymous (dS) substitutions of the PCGs were significantly larger than zero and smaller than 1, suggesting that all genes were experiencing purifying selection. The Complex V gene (atp6), the Complex III gene (cob) and the Complex IV genes (coxI-III) were under strongest selection, while the Complex I genes (nad1-6) were under weaker selection. These differences in strength of selection are reflected in the near horizontal relationship between dN and dS for the genes under strong selection (atp6, cob, coxI-III) and in a diagonal relationship for the genes under weaker selection (Fig. 4). The short genes from the Complex I genes (nad3 and nad6) exhibited a more scattered pattern. The Complex I genes had a higher synonymous substitution rate, with dS-values exceeding 0.5. Finally, the closest related species (PmI and PmIV) clearly accumulated less nonsynonymous mutations compared to the more distantly related species (Fig. 4).

Discussion
Free-living marine nematodes are characterized by a very high global species diversity, and nematode assemblages typically show high local diversity (Heip et al. 1985). Much work has focused on interactions among species and between species and their environment, and on the importance of nematodes for the functioning of marine benthic systems (Nascimento et al. 2012), yet only little is known on the evolutionary processes that have mediated this diversity. In addition, the presence of cryptic species may substantially complicate patterns of diversity. In marine nematodes, cryptic species have been detected based on diverging lineages in the mitochondrial COI gene and ribosomal nuclear genes (Derycke et al. 2005(Derycke et al. , 2007, and subsequent studies have demonstrated that these cryptic species show subtle morphological differentiation (Derycke et al. 2008a(Derycke et al. ,b, 2010, reproductive isolation , and ecological differences (De Meester et al. 2011, 2012aVan Campenhout et al. 2014). Our results now demonstrate high synteny in the mitochondrial genomes of the four species, but differences in gene length, in number of intergenic regions, and in the use of start/ stop codons between some species.
High synteny in the mitochondria of the L. marina cryptic species complex A survey in 62 nematode mitochondrial genomes, of which 54 species belong to the same order as our Litoditis marina species, revealed 25 different gene arrangement patterns (Liu et al. 2013). Variability in gene order is often due to a different location of the tRNA genes (Gissi et al. 2008), but this was not the case for the four Litoditis species. The gene order in the mitochondrial genomes of all four cryptic species was identical and corresponds to the typical gene order of several other rhabditid nematodes (Sun et al. 2014). The four mitochondrial genomes were highly AT-rich and codon usage was biased toward AT-rich codons, a pattern that has been frequently observed in Nematoda (Rota-Stabelli et al. 2010). The use of start and stop codons was highly comparable to that observed in other Rhabditina (Hu et al. 2002), including Caenorhabditis elegans (Okimoto et al. 1992), but different usage between the Litoditis group and C. elegans was observed for the nad4, cob, and nad1 gene. Furthermore, the use of incomplete stop codons such as T and TA has frequently been observed in Nematoda (Okimoto et al. 1992;Hu et al. 2002;Liu et al. 2013), but this was not the case for the Litoditis species. The use of the ATG start codon by PmIII for the initiation of the nad5 gene is unusual for metazoan mitochondria, but has also been reported in Ancylostoma duodenale (Hu et al. 2002), also a member of the Rhabditina. Transfer RNA in nematodes very often lacks one or two arms and shows replacement loops instead (He et al. 2005;J€ uhling et al. 2012). This was also observed for the L. marina species complex, where twenty tRNAs had a replacement T stem-loop, and the two tRNA Ser had a replacement D stem-loop. The annotated tRNAs further showed an overlap with the adjacent genes. Although overlaps up to 14 nucleotides have been recorded in enoplean nematodes (J€ uhling et al. 2012), the largest overlap in Litoditis was tRNA Thr with 30 nucleotides. Overlapping tRNAs can cause problems in tRNA processing, because they have a specific tRNA precursor and the upstream tRNA can miss the overlapping nucleotides (Reichert et al. 1998; M€ orl and Marchfelder 2001). The ribosomal RNAs of the Litoditis species showed only minor differences in length, and was comparable to the length of the ribosomal RNAs of C. elegans (697 and 953 bp, respectively) and of other nematodes (Hu et al. 2002).
The mitochondrial genomes of all four species lack the atp8 gene, which encodes for a core subunit of the F0 domain of the ATPase (da Fonseca et al. 2008). The absence of the atp8 gene has been observed in other nematodes as well (Okimoto et al. 1992;Montiel et al. 2006;Sultana et al. 2013;Leung et al. 2013), while a putative form has been found in the nematode Trichinella spiralis (Lavrov and Brown 2001). Genes encoding for the ATPase complex show a tendency to be lost in the fastevolving metazoan mtDNA (Gissi et al. 2008). Despite the high synteny in the mitochondrial genomes of the four cryptic species, our mitogenomic approach also revealed differences in genome length, gene length, nucleotide overlap between adjacent genes, and start and stop codon usage mostly between PmIII and the three other species.
The mitochondrial genome of PmIII shows several differences compared to PmI, PmII, and PmIV Differences in mitochondrial structural organization may have strong impacts on the physiology of the organism because mitochondria have an important role in energy metabolism. It was recently shown that a single nonsynonymous amino acid change in the mitochondrial cytochrome oxidase c subunit 1 (coxI) gene of geographical isolates of C. elegans substantially affects their longevity when exposed to different temperatures, which suggests that mitochondrial energy metabolism may be critical to adapt to environmental changes related to temperature (Dingley et al. 2014). Although we lack information on the longevity of these four species, recent experimental work shows that PmIII has a higher instantaneous fecundity than the other three species and performs better at higher temperature (De Meester et al. 2015). The coxI gene of the four species shows four nonsynonymous amino acid changes, three of which occur between PmIII and the three other species. PmIII further contains a substantially longer nad3 gene, which is part of the electron transport chain in Complex I (da Fonseca et al. 2008). Depending on their location and size, insertions and deletions can significantly alter the structure and function of proteins, and may be related to adaptation to particular environmental conditions (Wang et al. 2009). It is therefore possible that the structural differences in mitochondria may, at least in part, explain the biological differences between some cryptic species.

Phylogeny of the PCGs supports PmI and
PmIV as sister taxa and shows a Miocene origin for the Litoditis marina species complex The phylogenetic relationships obtained with the PCGs are consistent with previous published phylogenies of the Litoditis species complex based on nuclear data (Derycke et al. 2005(Derycke et al. , 2008a and recovered PmI and PmIV as sister taxa. Although relationships based only on the coxI gene have put PmII forward as the earliest diverging species (Derycke et al. 2005), our mitogenomic approach now supports PmIII as the earliest diverged species. This is in agreement with the large number of structural differences observed in the mitochondria of PmIII and with the high similarity in mitochondrial structure of PmI and PmIV. The presence of Wolbachia could affect the mutation rate and fixation of mutations due to sweeps in the mtDNA (e.g., Shoemaker et al. 2004) and can in this way affect phylogenies. In nematodes, Wolbachia is commonly found in filarial nematodes, but they have not been observed in secernentean nematodes (Bordenstein et al. 2003). Our data further support the lack of Wolbachia in Litoditis, a secernentean genus.
The most recent common ancestor for the Litoditis marina species complex was estimated at 16 MYA during the Miocene. PmII diverged 10.5 MYA from PmI and PmIV, while the latter two diverged 6.5 MYA. These dates are much older than the last glaciations and tectonic activities in the Atlantic Ocean with the formation of the North Atlantic-Arctic gateway and the Mediterranean Sea (Stoker et al. 2002;Harzhauser and Piller 2007), but they cannot be linked to mass extinction or other well-known geographical events. These old speciation events clearly show that cryptic species are not necessarily of recent origin. It is therefore unlikely that their similar morphology is caused by a lack of sufficient time to incorporate morphological differences. Morphological changes among species are low when strong selection on behavioral or physiological characteristics for adaptation to a specific host (Schonrogge et al. 2002) or environment are required. The Litoditis species thrive on decomposing algae in the intertidal environmental, where they are exposed to strong fluctuations in temperature and salinity. The ephemeral nature of the algal habitat requires the ability to rapidly colonize and reproduce to quickly establish viable populations. Litoditis is thus very likely to be subjected to strong environmental pressures. Visual recognition between species is often hampered in the marine environment and more pressure toward physiological distinctions (e.g., chemical characteristics) is likely to be important. We found evidence for purifying selection on all mitochondrial PCGs. Purifying selection is necessary for mitochondrial genes to maintain function (Rand 2001), although there are some reports of positive selection as well (summarized in (Castellana et al. 2011)). The coxI-III, cob, and atp6 genes are highly conserved in all four species, with strong purifying selection. The more closely related the species are, the higher the conservation in the genes between the species. The nad1-6 genes show higher variation with dS-values exceeding 0.5 in the genes with short sequences (nad3 and nad6). These genes have a higher rate of mutational saturation, which may lead to inaccurate dN/dS ratios. An earlier study in which dN/dS ratios were estimated from 347 complete vertebrate mt genomes showed that purifying selection was strongest for genes that encode subunits with crucial functions in the respiratory chain (RC), for example, cob and the coxI-III genes (Castellana et al. 2011). It has been put forward that the accumulation of non-neutral mitochondrial genetic variation within populations might play a role in speciation, with negative selection on mito-nuclear interactions (Dowling et al. 2008). These mito-nuclear interactions result in coevolution, with mtDNA mutations acting as drivers of adaptations in the nuclear genome (Dowling et al. 2008).

Conclusion
Despite millions of years of evolution, the mitochondrial genomes of the four cryptic Litoditis marina species show no variation in gene order or in gene content. However, many synonymous mutations have occurred, and structural differences between some species were present. Speciation of the cryptic species was estimated to have occurred in the Miocene, illustrating that the morphological similarity is not caused by insufficient time to evolve. Instead, differences in mitochondrial genome structures have accumulated between the earliest diverged species (PmIII) and the other three species, which may be linked to the different environmental adaptations observed in these species.