Hyper-Cryptic Marine Meiofauna: Species Complexes in Nemertodermatida

Nemertodermatida are microscopically small, benthic marine worms. Specimens of two nominal species, Sterreria psammicola and Nemertinoides elongatus from 33 locations worldwide were sequenced for three molecular markers. Species delimitation and validation was done using gene trees, haplotype networks and multilocus Bayesian analysis. We found 20 supported species of which nine: Nemertinoides glandulosum n.sp., N. wolfgangi n.sp., Sterreria boucheti n.sp., S. lundini n.sp., S. martindalei n.sp., S. monolithes n.sp., S. papuensis n.sp., S. variabilis n.sp. and S. ylvae n.sp., are described including nucleotide-based diagnoses. The distribution patterns indicate transoceanic dispersal in some of the species. Sympatric species were found in many cases. The high level of cryptic diversity in this meiofauna group implies that marine diversity may be higher than previously estimated.


Introduction
More than 70% of the earth's surface is covered by oceans, and sediment covers most of the ocean floor. Marine infauna thus inhabits one of the earth's largest ecosystems. Sediment meiofauna is a diverse assemblage with representatives from many animal phyla. Despite the vast size of the marine benthic ecosystem, the marine meiofauna is poorly known, and even in well-studied areas numerous undescribed species exist [1][2][3][4]. Nominal meiofauna species are often reported to have cosmopolitan distributions in concordance with the ''Everything is Everywhere (EiE)'' hypothesis stating that animals below 1 mm body size are easily dispersed. EiE was originally applied to microorganisms [5] and later extended to organisms up to 1 mm size [6,7]. However, species identification of meiofauna requires time-consuming microscope studies, which is often only possible when a specialist brings equipment to the field to examine live specimens. Such detailed taxonomic studies have shown a high level of endemicity for some groups, e.g. Platyhelminthes and Acoela [3], thus contradicting the EiE hypothesis; whereas other groups such as gastrotrichs of the genus Turbanella seem to conform to a pattern of large distributions [8].
The diversity of the marine worms of the taxon Nemertodermatida that are part of the meiofauna in clean sandy sediments was reviewed by Sterrer [9] who recognized eight broadly circumscribed species with a potential for further subdivision as some of them were known only from few specimens, many of which were incomplete. Morphological identification of nemertodermatid species is complicated by the fact that a large number of specimens are juveniles where the diagnostic reproductive organs cannot be studied. In total Sterrer [9] reported that 229 specimens of nemertodermatids were studied by him since 1964. The nominal species with the largest distribution range was Sterreria psammicola Sterrer, 1970. Sterrer studied 43 specimens of S. psammicola from the North Sea area, the Mediterranean, Caribbean, Australia and Papua New Guinea and considered it ''remarkably homogeneous throughout its global distribution range'' and regarded Nemertoderma rubra (Faubel 1976) [10] as its junior synonym. There is, however, some morphological variation in this cosmopolitan species, most apparent in the pigmentation, which can range from non-existent, with the worms appearing glossy silvery, over a narrow, often only anterior, reddish or brownish ''spinal stripe'' to a more or less uniform bright red colour (Fig. 1a, h). Pigmented and unpigmented specimens have been recorded from the same site, e.g. around the island of Helgoland, North Sea. The nominal species Nemertinoides elongatus Riser, 1987 [11], which is known only from relatively few specimens, is similar in shape to S. psammicola and juvenile specimens of the two species cannot be distinguished although adults differ in reproductive anatomy and the morphology and distribution of epidermal gland cells.
While current taxonomy suggests that Sterreria psammicola has a cosmopolitan distribution, there are biological factors that indicate limited dispersal ability and consequently a high degree of endemism in these small interstitial marine worms: none of the known nemertodermatid species have dormant eggs or a planktonic stage; small and fragile juveniles hatch from thinshelled eggs shortly after they have been deposited [12].
A recent estimate of marine eukaryote biodiversity based mainly on expert opinion concluded that there may be 0.7-1.0 million marine species including the 226 000 currently known nominal species. The proportion of cryptic species remaining to be identified was approximated to range between 11% and 43% of the currently known number [13]. Here we follow Bickford et al. [14] in regarding as cryptic those species that are or have been classified as the same nominal species due to morphological similarity. Some groups, including Nemertodermatida, were considered too poorly known to allow an estimate of the incidence of cryptic species by Appeltans et al. [13]. Costello et al. [15] estimated the number of marine species based on the rate of descriptions of new species using data from the World Register of Marine Species (WoRMS) and concluded that there are 300 000 marine species. The latter study did not discuss the effect that undetected cryptic species would have on the estimate. Neither of the two studies defined their concept of species, instead a ''legacy species concept'' based on the numbers entered into WoRMS was used.
Currently, developments in DNA-sequencing technology and bioinformatics are unleashing the potential for broader and deeper sampling of marine biodiversity. Poorly known meiofauna taxa that may exhibit low morphological complexity or be fragile may thus become better known through metagenetic studies such as those by Fonseca et al. [16]. Next-generation sequencing is also likely to reveal additional diversity in the form of cryptic species; see Emerson et al. [17] for an example from soil fauna. Metagenetic studies of diversity will be immensely more valuable when a populated database of sequences from known species with revised taxonomy is available.
Here we aim to test the hypothesis that the nominal species Sterreria psammicola and Nemertinoides elongatus are complexes of cryptic species. Our data is based on 172 specimens that were collected during seven years from 33 different locations (Tab. 1). We sequenced complete or near-complete ribosomal large and small subunit (LSU and SSU) genes and a fragment of the protein coding Histone 3 (H3) gene, and computed separate gene trees under Maximum Likelihood and Bayesian approaches as well as parsimony networks and pairwise distances to generate primary species hypotheses. Clades identified as putative species were tested for genetic isolation using a multilocus Bayesian approach with the software BP&P [18] to generate secondary species hypotheses. Clades with at least three specimens that were supported in at least two of the three gene trees, present as separate haplotype networks under statistical parsimony, had an averaged interspecific pairwise distance at least twice the averaged intraspecific distance, and that were validated by multilocus Bayesian analysis, are formally described and named in this paper. We operate with a species concept in accordance with the ''unified species concept'' of de Queiroz [19] emphasizing that species are independently evolving lineages that can be diagnosed in a multitude of ways.

Permits
Taxa used in this study are interstitial invertebrates, which do not need special sampling permits, as they are not subject to regulations of species protection and are collected within small amounts of sediment.
For sampling around Helgoland, Germany, at Waimanolo, Hawaii, in Norway, Sweden and most of the Mediterranean, no specific or additional sampling permits for the collection of small amounts of marine sediments were required. Geographic coordinates for each site are given in

Specimens
Specimens were extracted from sediments using isotonic magnesium chloride solution [20] and identified under a dissecting microscope sometimes in combination with a compound microscope. Specimens were photographed using a compound microscope, if possible equipped with differential interference contrast optics, before fixing in ethanol or RNAlater. Their microscopic size necessitates use of whole specimens for DNA extraction. To ensure a direct link between morphology and gene sequences all type specimens were photographed prior to preservation for DNA extraction and images are deposited as illustrations of the type material, see table 2 for museum and genbank accession numbers. For the description of the position of morphological characters, a relative scale (U) is used with the anterior tip of the animal corresponding to 0 U and the posterior tip to 100 U [21]. Measurements, however, are difficult to take as animals seldom lie straight and relaxed for a sufficiently long time and in many cases specimens are incomplete, as the worms are fragile.
DNA extraction, amplification and sequencing DNA was extracted using the Qiagen Micro Tissue Kit. The microscopic size and corresponding low yield of extracted DNA from the specimens as well as the unavailability of prior sequence data severely limited the choice of nucleotide markers. We were able to consistently amplify and sequence rRNA genes as well as the nuclear protein coding Histone 3 gene. The large ribosomal subunit gene was obtained from 168 specimens with an alignment length of 3583 bp, the small ribosomal subunit gene from 166 specimens (1792 bp) and H3 from 106 specimens (328 bp). All markers were amplified and sequenced using several different primer combinations (Tab. 3), and, in the case of SSU, a nested PCR approach.
Sequence editing, alignment (MAFFT [22]), translation into amino acids and checks for open reading frames were performed using the Geneious Pro 7.0.4. software package created by Biomatters available from http://www.geneious.com. The alignments were tested for random similarity with the program Aliscore [23,24] using the default settings. jModeltest v. 2.1.1. [25] analyses were performed for each dataset in order to test the datasets for the use of the proportion of invariable sites (I, propinvar) and the rate variation across sites (G) and to obtain values to set useful priors. Evolutionary neutrality of the coding gene H3 was tested using Tajima's D calculated with the software MEGA 5 [26]. Saturation of the H3 gene was detected through plotting the uncorrected pdistances versus the phylogenetic distance using an R-script [27].

Phylogenetic ''species discovery''
The Geneious package (v. 7.0.4.) was also used to calculate pairwise distances between sequences within and between putative species. For this the LSU and SSU alignments were trimmed by eye to 2009 bp and 1502 bp respectively in order to have sequences of similar lengths but keep most of the information. Those specimens represented by less than half of the alignment length were excluded (s. Supplementary table ST1 for details).
Parsimony haplotype networks were computed using the software TCS 1.21 [28] with the reduced and trimmed datasets for LSU (further reduced to 154 specimens and 2009 bp) and SSU. Gaps were considered a fifth state. For relatively fast evolving mitochondrial genes, a 95% threshold has been shown to recover known species reliably [29]. To account for a slower evolutionary rate the connection limit was set to 98% for the LSU and SSU genes; an additional analysis with the connection limit of 90% was performed for the higher resolving Histone 3 dataset. Maximum Likelihood (ML) and bootstrap support calculations were performed by raxmlGUI [30] using the GTR+G+I evolutionary model and the rapid bootstrap algorithm with 1000 bootstrap reiterations.
Bayesian analyses were performed using the program MrBayes 3.2.1. [31]. No evolutionary model was set and the program was allowed to sample the entire model space of the GTR model by defining nst = mixed. The proportion of invariable sites and G were applied with the prior set to shapepr = Uniform(0.05,1.00) for SSU and LSU and shapepr = Uniform(0.05,2.00) for H3; the pinvarpr was left at the default. Analyses were stopped when the standard deviation of split frequencies was between 0.01 and 0.05, indicating sufficient convergence and a relative burn-in of 25% was used. No concatenated analyses for all three genes combined were conducted. This would conceal incongruences between the gene trees and therefore possibly lead to subsequent errors in the validation of species using BP&P [32].

''Species'' validation
Those clades that consisted of at least three specimens, showed an averaged interspecific pairwise distance at least two times higher than the intraspecific averaged pairwise distance (relative threshold distance [34]), formed separate parsimony networks and were present in at least two of the three gene trees, were tested using a multilocus Bayesian approach with the program BP&P to generate secondary species hypotheses [18,35], species validation sensu [36]. The program relies on a user-defined tree and only tests for the presence of nodes in the input-tree; the input of an incorrect guide tree will corrupt the results [32]. In order to create unambiguous input trees the dataset was divided into three subsets and the putative species Sterreria martindalei n.sp. and Sterreria papuensis n.sp. were excluded (different colours in Fig. 2, 3) because the gene trees could not resolve all deeper nodes with high support. Both excluded species, however, are highly supported in all species discovery methods, thus we think that further validation in these cases was not necessary. The subgroups within Sterreria variabilis n.sp. were not validated because of the unresolved topology (polytomies) of the group. Two analyses with the gamma priors set to G(1, 100) and G(1, 1000) for the population size h and G(1, 100) and G(1, 1000) for the root age t were conducted while the other divergence time parameters are assigned the Dirichlet prior ( [18]: equation 2). An additional analysis with an older root age with the G of h (1, 100) and the G of t (1, 10000) was also conducted.
The species we describe are diagnosed based on unique differences in the nucleotide sequences following Jörger and Schrödl [4] in addition to morphological diagnostic characters, which are provided where available.

Nomenclatural acts
The electronic edition of this article conforms to the requirements of the amended International Code of Zoological Nomenclature, and hence the new names contained herein are available under that Code from the electronic edition of this article. This published work and the nomenclatural acts it contains have been registered in ZooBank, the online registration system for the ICZN. The ZooBank LSIDs (Life Science Identifiers) can be resolved and the associated information viewed through any standard web browser by appending the LSID to the prefix ''http://zoobank.org/''. The LSID for this publication is: urn:lsid:zoobank.org:pub:A306F670-B4B4-4376-A859-48A9735E1593. LSIDs for new species are given

Results
When testing for random similarity between sequences, Aliscore highlighted 419 of 3583 aligned sites in the LSU dataset and 110 of 1792 sites in the SSU dataset. No random similarities were indicated in the H3 dataset. Consensus and best trees resulting from analyses of the original and Aliscore-filtered alignments had identical topologies with the exception of four specimens of Sterreria rubra that grouped with the specimen S7 in the Aliscore pruned LSU analysis. There were small differences in branch support when comparing original and filtered alignments (greatest difference 6% BS in the LSU dataset and 55% in the SSU dataset, but generally no or less than 10% BS difference). The gene trees shown as supplementary data are based on the original alignments ( Figures S2-7). A summary of the results in terms of putative species is given in figure 2, which shows a 75% majority rule consensus tree (MF75) of the LSU ML analysis performed with RAxML. Support (bootstrap for ML analyses, posterior probabilities for Bayesian analyses) for the 20 putative species (excluding outgroups) is shown to the right of the node (or in case of a few long branches over those) in the order LSU/SSU/H3 (for all gene trees see Figures  S2-7). The colours refer to the subsets used for species validation (green: Nemertinoides-group, red: European Sterreria-group, blue: extra-European Sterreria-group, orange: untested Sterreria species).
In the uncorrected pairwise distances matrix several groups with at least twice the intraspecific distance to their sister group could be identified (Tables S2-4). The LSU and SSU gene datasets each had 19 distinct putative species groups (excluding outgroups), and in the Histone 3 gene dataset we found 28 such groups. The intraspecific distances never exceeded 0.8% and 0.5% respectively in the LSU and SSU data partitions, with the exception of N4 with 1.2% in the LSU dataset. In the LSU dataset the single specimen representing putative species S7 was excluded from the pairwise distance analysis because the sequence was too short. In the SSU dataset the species Sterreria ylvae n.sp. and S. monolithes n.sp. could not be distinguished from each other (averaged interspecific pairwise distance 0.2%). In the H3 partition intraspecific distances reached 9.9% in Sterreria variabilis n.sp. Of the 27 groups in the H3 data, 15 correspond to the same putative species as seen in the LSU and SSU gene datasets (Table 4). Eight groups in the H3 dataset did not correspond to putative species supported by the other two genes. This may be a saturation artefact ( Figure S1).
The TCS software defined different numbers of parsimony haplotype networks for each of the three loci. Analyses of the LSU, SSU and H3 gene datasets with a connection limit of 98% found 25, 20 and 42 networks respectively (excluding outgroups, Fig. 3a, b). When the Histone 3 gene analysis was relaxed with a connection limit of 90% only 30 networks were found (Fig. 3c). In the LSU dataset N. glandulosum n.sp., S. papuensis n.sp., S. psammicola and S. variabilis were recovered as two and three separate networks respectively. Putative species S7 was excluded from the dataset due to its short sequence. In the SSU analysis, S. boucheti n.sp. and S. ylvae n.sp. were recovered as one network with two steps between the two species. One specimen of the diverse S. variabilis n.sp. formed a separate network not connected to the other specimens of the species. The H3 gene analyses split S. lundini n.sp. and S. papuensis n.sp. into two networks each and N. elongatus into three different networks. S. rubra was recovered in seven networks most of them consisting of only one or two specimens, corresponding with the observed pairwise distances. S. variabilis n.sp. formed five networks and one network connecting with S. boucheti n.sp. S. ylvae n.sp. and S. monolithes n.sp. formed one network connected by ten steps. In summary the network assemblages discovered with TCS are highly congruent with the groups identified in the pairwise distance matrix between genes, especially in the LSU and SSU genes. Tajima's D for the H3 dataset is D = 1.931985, which indicates that this marker underwent neutral evolution. The H3 saturation test indicates saturation (S1).
The gene trees of the three loci estimated with RAxMl (best tree with bootstrap values) and MrBayes are not resolved in the deeper topology but consistently support the same putative species (Fig. 2,  Figures S2-7). In the LSU and SSU genes 20 putative species (excluding outgroups) can be identified and in the Histone 3 dataset 27 such groups are supported. The groups identified in the gene trees are identical or highly congruent with the groups identified in the pairwise distance analyses and by the haplotype networks.
In the H3 gene trees thirteen of the 101 ingroup specimens are recovered as one clade splitting from a basal trichotomy. These specimens belong to S. lundini n.sp., S. papuensis n.sp. and S. rubra. This grouping can be interpreted as a saturation artefact. Exclusion of these specimens from the analyses did not change the composition of other tip groups (putative species). The recovered putative species in the gene trees, other than S. lundini n.sp. and S. papuensis n.sp., are consistent with those identified from the pairwise distances and haplotype networks. Table 4 summarizes the support for identified groups over all methods of discovery.

Species validation
The BP&P analyses with less informative priors or an older root age supported all eleven putative species tested (Fig. 4). The analysis with an informative prior supported all putative species except Sterreria lundini n.sp., S. psammicola, S. rubra and S2 (Fig. 4).
All 35 specimens with any kind of rosy, brownish to bright red colouration recorded were found in the species Sterreria rubra and clade S2 together with four uncoloured specimens and seven specimens for which no colour data was recorded.
Seven out of the eight Sterreria species are not globally distributed, with two species being limited to Hawaii and Papua New Guinea respectively (Fig. 5). Only one species, S. variabilis n.sp., includes specimens from Hawaii, Papua New Guinea, New Caledonia, Bermuda, Portugal and the Mediterranean.
All three Nemertinoides species are distributed along the European coastline from the Mediterranean, via Portugal and the North Sea to the Swedish West Coast (Fig. 5). Within each clade however, no patterning corresponding to geographical distances could be observed. The same is true for the putative species indicated with red ( Figs. 2 and 5), which also consist of European specimens and show no geographic pattern.
Clades meeting the above mentioned criteria and were supported in the BP&P analyses (except S. martindalei n.sp. and S. variabilis n.sp.) are formally described, with the exception of clade S2; the specimens of this clade lack photographs or sketches, prohibiting formal description of the species (Tab. 4).

Discussion
Our dataset multiplies the available nucleotide sequence data for Nemertodermatida more than 100 times. The analyses of nucleotide sequence data from the three genes for LSU rRNA, SSU rRNA and Histone 3 identified twelve well supported species among the collected specimens from the two nominal species Nemertinoides elongatus and Sterreria psammicola (Tab. 4). N. wolfgangi n.sp. and S. lundini n.sp. are supported in all genes and analyses but are paraphyletic in the SSU gene trees. The unified species concept defines species as a ''separately evolving lineage segment'' [19]. SSU has been shown to underestimate species diversity in meiofauna [37]. We therefore conclude that both species are independently evolving lineages warranting formal description as species with incomplete lineage sorting in the SSU gene.
Twelve described species, however, clearly is an underestimation of the true biodiversity in this group of Nemertodermatida (Fig. 2). Even in the material from the Mediterranean, the most densely sampled geographical area, there are several clusters of less than three specimens consistently grouping together in all analyses. These clades represent additional cryptic diversity but we refrain from formally naming such poorly sampled putative species here. This taxonomic undersampling is of course even more drastic outside the Mediterranean. A different form of undersampling emanates from the limited dataset that we have acquired: inclusion of additional molecular markers would have boosted the potential to detect additional cryptic species. Our conservative approach, based on three molecular markers, still raises the number of named species of Nemertinoides and Sterreria from two to twelve.
With the material available to us, morphological distinguishing characters could not be identified a priori for all of the herein described species while studying live specimens, but it is possible that such features will be discovered a posteriori if more material becomes available. We suspect that this is a matter of the level of detail in the morphological investigation, which could be extended from light microscopy to CLSM or electron microscopy in search of additional characters. Even if some of the new species remain diagnosable only based on nucleotide sequences, it is important to recognize such cryptic species in order to appreciate biodiversity, plan management of conservation, and understand ecosystem function [14] (with references). This is especially true if much of species diversity is constituted by cryptic taxa, as is evidently the case in Sterreria and Nemertinoides. It has been argued that species delimitations based solely on nucleotide sequence data would lead to taxonomic instability and confusion as well as taxonomic inflation [38]. However, multilocus coalescent-based methods for species delimitation are firmly grounded in evolutionary theory and population biology, and since these methods are based on explicit probabilities they can be considered more objective than traditional character based taxonomy and allow greater comparability between species [39,40]. Furthermore, if nucleotide-based species diagnoses were implemented, juvenile specimens as well as fragments of specimens, a large proportion of the nemertodermatid specimens encountered, would be available for the study of the diversity of this group.
Adams et al. 2014 [41] defined hyper-cryptic species as nominal species that actually consist of four or more valid species. Our application of molecular species discovery tools have revealed that the two nominal species Nemertinoides elongatus and Sterreria psammicola are hyper-cryptic as they are composed of at least 20 separate species-level clades. Nine of these will be formally described and named below, thereby doubling the number of nominal species of Nemertodermatida. There is no reason to believe that Nemertodermatida are unique in their extensive cryptic diversity: analyses of nucleotide sequence data have unravelled cryptic and hyper-cryptic species within many other groups of marine invertebrates. A case in point is the nominal polychaete species Eumida sanguinea (Ö rsted, 1843) which was studied by Nygren and Pleijel [42] who identified eight cryptic species among specimens assigned to E. sanguinea and named seven of them using nucleotide-based diagnoses. There are a number of additional cases of cryptic diversity in other polychaete taxa, e.g. the ''cosmopolitan'' fireworm Eurythoe complanata (Pallas 1766), which was found to consist of three species [43], and Notophyllum foliosum (Sars 1835), which was found to consist of two species [44]. There are relatively few studies at this level of taxonomic resolution in marine meiofauna, but cryptic species have been identified in the flatworm genera Pseudomonocelis [45,46] and Monocelis [47]. A noteworthy example is the ''cosmopolitan'' flatworm Gyratrix hermaphroditus Ehrenberg, 1831, where studies of karyotype and fine morphology revealed eight separate species in Australia [48], two separate species in the North Sea and the Mediterranean and two separate species at the French Atlantic coast [49]. Leasi et al. [50] used the coalescentbased GMYC algorithm [51] to analyse Cytochrome oxidase subunit I sequences from specimens of the rotifer Testudinella clypeata (Müller, 1786) and found seven cryptic species. Our results further corroborate the hypothesis of the oceans as a hotspot for cryptic diversity put forward by Bickford et al. [14] and exemplified above. Cryptic species occur in all animal groups and they are being identified and described at an accelerating rate Table 4. Summary of the results of the different species identification methods (per genetic marker) and the multi-locus species validation (BP&P). The number of specimens is given per gene in the order LSU, SSU and H3. The uncorrected percentage of pairwise distances was calculated on the same reduced dataset as was used for the calculation of the parsimony networks; the lowest average percentage of pairwise distance between any two clades within the datasets versus the averaged intraspecific pairwise distance is shown. The parsimony networks were calculated with a threshold of 98% for LSU and SSU and 90% for H3; numbers specify the steps between the two clades in case they were recovered in one network, an asterisk indicates that a given clade was recovered in more than one network. Gene trees were calculated using RAxML and MrBayes, support is given in this order; non-monophyletic clades are indicated with a hash in superscript ( # ). The program BP&P relies on a true guide tree, which due to the low taxon sampling we could not provide. We therefore decided to test only species with at least three specimens, thus some of the clades in the dataset were not tested and will not be named in this paper (except for S. martindalei and S. papuensis, which were well supported and described without BP&P validation). The posterior probabilities of the three different BP&P analyses are shown in the table. Due to the high variability of S. variabilis clade in all analyses, the clusters within this clade were tested for all analyses with the H3 gene. However, as this part of the tree was not resolved in either phylogenetic analysis, so no BP&B validation could be performed. A dash indicates missing data, or not tested clades. doi:10.1371/journal.pone.0107688.t004 [14,52,53]. We hypothesize that the amount of cryptic diversity in meiofauna is far higher than the 11%-43% estimate proposed by Appeltans et al. [13]. Consequently, the approximation of total marine diversity at 0.7-1.0 million species is likely to be an underestimate.

Species/clade
As noted above, identification and taxonomic study of Nemertodermatida requires specialized methods and access to live specimens, which may be part of the explanation for our fragmentary knowledge of their diversity. This is true also for other groups of meiofauna, especially fragile groups such as Acoela, Platyhelminthes and Gastrotricha that cannot be easily preserved for future identification. Application of metagenetic methods e.g. [16,54,55], where DNA is extracted from sediment samples followed by PCR amplification of a selected marker, pyrosequencing and bioinformatic analysis, has potential to change this as the morphological identification stage is eliminated. This procedure is considered cost-effective in comparison to traditional methods where a number of taxonomists would have to study each sample [56]. It will also provide a more complete snapshot of diversity, as juvenile or damaged specimens would contribute to the results. Identification of species through the metagenetic approach requires a populated database of reference sequences from specimens that were identified by a specialist, such as those provided in this study.

Habitat and Biogeography
All our specimens were found in depths between 1.5 m and 37 m in sand that reached from coarse to very fine with low to moderate organic content. However, we did not observe any differences in habitat specific to the identified species; if present they are clearly subtle. The nemertodermatid species Nemertoderma westbladi and Meara stichopi, outgroup species in this study, occur in mud down to depths of 600 m. Sampling of nemertodermatids from deeper sediments has only been done in very few locations, mainly in the North-East Atlantic. Thus the existence of deep-water species of Sterreria and Nemertinoides cannot be ruled out.
No fossils identified as nemertodermatids are known, but according to both of the two currently competing hypotheses on the phylogenetic position of Nemertodermatida [57,58], the group is as old as the Cambrian explosion, or even predating it. When trying to explain the current distribution of such an old clade as the Nemertodermatida, and taking into account their biology with the poor capacity for dispersal that it implies, a first hypothesis may be to invoke vicariance, explaining the patterns by continental drift in combination with speciation. However, the vicariance hypothesis cannot explain the presence of littoral Nemertodermatida on younger and isolated islands such as Hawaii and Bermuda. O'ahu island, where our Hawaiian specimens were collected, is three million years old [59]. It should also be noted that the high estimated age of Nemertodermatida, deduced from their phylogenetic position, pertains to the clade as a whole. The age of the nemertodermatid crown group, which includes the recent species of Sterreria and Nemertinoides, cannot be determined in the absence of any calibration points, but it is likely to be much younger. Clearly, dispersal is the only feasible explanation for the presence of Sterreria on O'ahu and Bermuda. The possibility remains that isolated young islands were colonized by deepwater nemertodermatid species, although currently available evidence seems to favour dispersal from distant shallow habitats, as no deep-water Sterreria specimens are known. The phylogenetic hypothesis derived from Bayesian analyses of the ribosomal datasets indicates that O'ahu was colonized at least twice as Sterreria ylvae n.sp., Sterreria martindalei n.sp. and Sterreria variabilis n.sp. are not each others closest relatives. Current Table 5. Molecular Diagnostic character of all newly described species in the three genes used in this study.

LSU SSU H3
sampling of Sterreria on O'ahu was restricted to one site. Extended sampling of nemertodermatids in the Hawaiian archipelago, where the islands are of different ages ranging from 28 Myears to 400 000 years [59], would allow an estimate of rates of colonization and speciation within Sterreria as well as indicating the relative importance of dispersal and speciation in nemertodermatid diversity. Similar studies in Macaronesia, with its volcanic islands of different ages and degrees of geographical isolation, would also shed light on the genetic distinctness of the Bermudian species and the transatlantic dispersal.
Our results show that Nemertodermatida mostly do not conform with the EiE hypothesis. The supposedly wide-ranging Sterreria psammicola and Nemertinoides elongatus both consist of complexes of cryptic species. Some of the species, e.g. Nemertinoides wolfgangi n.sp., show a distribution pattern restricted to one ocean, in this case the Mediterranean, which is, as noted above, the best represented area in this study. However, other species have more extensive distribution areas, such as Nemertinoides elongatus and Sterreria rubra, ranging from the Adriatic, through the Mediterranean via Portugal and the North Sea into the Skagerrak without exhibiting any apparent genetic structure. This indicates considerable dispersal ability in these interstitial worms. Outside Europe, the findings in the Madang lagoon, Papua New Guinea, are particularly striking: the 25 animals collected from four adjacent localities (less than 10 km distance) belong to six different species (S. papuensis n.sp., S. variabilis n.sp., S. boucheti n.sp., S. monolithes n.sp., S7 and P3). Of those, S. papuensis n.sp. and S7 are more closely related to species with European distribution, than to species from geographically closer localities. This is clearly not consistent with an isolation by distance pattern and again indicates dispersal.
Only one species in our dataset appears to be truly cosmopolitan: Sterreria variabilis n.sp. However, since the Histone 3 gene splits this nominal species into geographically structured clades, the existence of yet another unresolved species complex is possible. Remarks: The far anterior position of the male pore in N. elongatus with posterior ovaries were the main arguments for the naming of a new genus by Riser [11]. However, the position of the male pore in relation to body length proved to be not informative on genus level, as opposed to the relative position of the ovaries reaching posterior of the male opening.  (7), and the Adriatic (3). More detailed information about individuals and further photographs are accessible at http://acoela.myspecies.info/, the scratchpads database for Acoela and Nemertodermatida.

Taxonomic part
Description: Up to 6 mm long, region anterior of statocyst slimmer than rest of body. Statocyst at U4. Large bottle glands in epidermis in anterior half of body. Mouth at U25. Male opening at U40 (Fig. 6e); false seminal vesicle (fsv) directly anterior to this (Fig. 6d). Paired ovaries extend from posterior of fsv to posterior tip. Found in slightly coarser sand from the intertidal to 30 m depth.
Remarks: This species conforms to the description of N. elongatus in Riser [11] with the position of the male opening at about U40 and the distribution of epidermal glands only in the anterior half of the body length. Species supported by all three genes in this study. The type material for N. elongatus was collected in the Western Atlantic (holotype at the Massachusetts coast and paratypes at the New Brunswick coast). Subsequent specimens identified as N. elongatus were collected in the Adriatic Results given as nodal support for all Nemertinoides species (green in Fig. 2), mainly European Sterreria species (red in Fig. 2) and the extra-European Sterreria species (blue in Fig 2). Support values are Bayesian posterior probabilities for the different analyses in the order G(1/100), G(1/1000) and old root age (G(1/100) and G(1/10000)). The dataset was split in order to avoid artefacts due to unresolved topologies in the gene trees and increase confidence in the input topologies. Only clades represented by more than two specimens were tested in order to increase confidence. doi:10.1371/journal.pone.0107688.g004 by Sterrer [9]. There is some uncertainty attached to the identification of this species, as we did not have access to specimens from the type localities; it is possible that our specimens represent a species different from the original N. elongatus.
Distribution: Western Atlantic, Swedish West coast, North Sea, southern Portugal, Mediterranean.
-Nemertinoides wolfgangi n.sp. Material examined: 9 living specimens in squeeze preparation collected mostly in summer between the years 2008 and 2013 in the French Mediterranean (1), the Tyrrhenian Sea (5), Sicily (2), and the Adriatic (1). More detailed information about individuals and further photographs are accessible at http://acoela.myspecies.info/, the scratchpads database for Acoela and Nemertodermatida.
Remarks: Species supported by all three genes in this study. Distribution: Mediterranean, southern Portugal, North Sea. Etymology: Glandula = latin for gland. This species, in contrast to the prior described N. elongatus, has epidermal glands also in posterior.
Remarks: The species Sterreria psammicola has been described originally within the genus Nemertoderma by Sterrer in Riedl (1970). Due to differences in epidermal structure and differences in the position of the reproductive glands Lundin [60] created a new monotypic genus Sterreria and placed Sterreria psammicola in this. The mouth has been observed only in few specimens, it is hypothesized to be a temporary feature as in Nemertoderma westbladi [61][62][63]. ''Male maturity seems to precede female maturity'' [9].
Remarks: Sterreria psammicola has been described from Croatia as a filiform worm living in shallow sandy sediments under the name Nemertoderma psammicola by Sterrer in Riedl [64]. The species was described based on some specimens in which ''a salmon-red longitudinal stripe is usually (italics by the present authors) present in the first third of the body length'' [64]. Specimens from this clade conform to the description without colouration and two specimens were collected near the type locality of Sterreria psammicola in Croatia. In the formal description in 1970 no type material was deposited. Species supported by all three genes in this study.
Description: Usually pigmented, rose, bright red or brownish; only in anterior part or all over body (Fig. 1a, b; 7 g-j). Frontal glands prominent reaching far behind statocyst; opening fanning out at the anterior tip; secretions rod shaped (Fig. 7i). Epidermal glands small, distributed densely especially in anterior third of body, but never as prominent as in Nemertinoides-species. Body surface appears ''scaly'' due to visible epidermal cell borders (Fig. 7j). Testes lateral. Male pore at U90; fsv just anterior to that (Fig. 7g). Ovaries paired, usually two mature eggs and several small oocytes (Fig. 7h). Posterior tip wide.
Remarks: In 1976 Faubel described the species Nemertoderma rubra from the islands Rømø and Sylt in the North Sea based on three specimens and stated that ''in transmitted light the species is coloured reddish'' [10]. In a revision of the taxon Nemertoder-matida, Sterrer [9] remarked that Nemertoderma rubra and N. psammicola are very similar and regarded them as one species, making N. rubra a junior synonym of N. psammicola. Due to the specimens of this clade conforming to the comprehensive formal description of Nemertoderma rubra and the clear statement that ''the species is coloured reddish'' we reinstate the junior synonym Sterreria (Nemertoderma) rubra. Ovaries are positioned further towards the posterior than described by Sterrer [9] and Faubel [10]. Species supported by all three genes in this study, but polyphyletic in the H3 gene tree (saturation artefact).
Distribution: North Sea, Swedish West coast, southern Portugal, Mediterranean.
-Sterreria lundini n.sp. Material examined: 12 living specimens in squeeze preparation collected mostly in summer between the years 2008 and 2013 in the North Sea (1), the French Mediterranean (1), the Tyrrhenian Sea (9), and Sicily (1). More detailed information about individuals and further photographs are accessible at http://acoela.myspecies.info/, the scratchpads database for Acoela and Nemertodermatida.
Typematerial: Holotype SMNH type-8634 (collection code UJ-08117). Male mature specimen collected near Castiglione della Pescaia, Italy, by Marco Curini Galletti 19 May 2008. Photographs of the holotype deposited at the Swedish Museum of Natural History, Stockholm.
Remarks: This species is paraphyletic in the SSU gene tree and polyphyletic in the H3 gene tree (saturation artefact).
Distribution: Mediterranean. Etymology: After Kennet Lundin, the first researcher formulating a comprehensive phylogenetic hypothesis for Nemertodermatida.
-Sterreria papuensis n.sp. Material examined: 10 living specimens in squeeze preparation collected in November 2012 near four different islands (Siar, Tab, Panab, Wanad) in Madang Lagoon, Papua New Guinea. More detailed information about individuals and further photographs are accessible at http:// acoela.myspecies.info/, the scratchpads database for Acoela and Nemertodermatida.
Typematerial: Holotype SMNH type-8637 (collection code PNG77). Immature specimen collected near Wanad island in November 2012 by Inga Meyer-Wachsmuth. Photos of the holotype deposited at the Swedish Museum of Natural History, Stockholm.
Description: Up to 1 cm long. Opaque, dirty rose. Frontal glands reach far towards posterior branching tree-like, opening fanning out. Mouth at U35. Male opening at U85. Adhesive structure in posterior. Borders between epidermal cells visible (''scaly'') similar to Sterreria rubra. Epidermal glands few but regularly distributed.
Remarks: Species supported by all genes in this study but is polyphyletic in the H3 gene tree (saturation artefact).
Etymology: papuensis = coming from Papua (New Guinea). This is remarkable as this species is closely related to otherwise exclusively European (putative) species.
Remarks: All three specimens of species had only one statolith in the statocyst. Aberrant numbers of statoliths can be observed frequently in other species of Nemertodermatida without taxonomic implications. However, it is unusual to find only specimens with an abnormal number of statoliths. All studied specimens were immature.
Distribution: Madang lagoon, Papua New Guinea. Etymology: monos = greek for single, lithos = greek for stone, referring to the single lithocyte in the statocyst in all three studied specimens.
Diagnosis: Statoliths appear to be encapsulated into separate membrane inside statocyst. Molecular character diagnosis in Description: Statocyst at U5. Overall larger, sturdier than Sterreria ylvae n.sp. Posterior two thirds opaque due to large and dense epidermal glands dorsally or ventrally, not laterally. Frontal gland opening central, secretions ellipsoid, reaching posterior till U35; clear border between frontal and epidermal glands. Ovaries between U50 and U75 (Fig. 8a).
Remarks: Species supported by all three genes in this study. Distribution: O'ahu Island, Hawaii, USA. and adhesive area (ad). e) Overview over male mature animal with large epidermal glands distributed more or less evenly over the whole body and mco in the posterior. f) Anterior part of an animal with epidermal glands, statocyst and central frontal gland opening (fp). g-j: Sterreria rubra. g) Posterior with two mature eggs (e) and mco. h) Overview over the same male and female mature specimen with mature eggs (e) and oocytes (oo). i) Anterior with statocyst and rod-shaped frontal gland secretions. j) Detail of the epidermis in the anterior of an animal with cell borders clearly visible, giving the animal a ''scaly'' appearance. k: Stererria papuensis n.sp. Anterior part of the holotype. All photographs are taken of live specimens in squeeze preparation. a, b, c, f and k are photographs of the holotype and neotype specimens respectively. The scale bars indicate 100 mm for each photograph. doi:10.1371/journal.pone.0107688.g007 Etymology: After Prof. Mark Martindale, who hosted U. J. in Hawaii and made sampling possible.
Typematerial: Holotype SMNH type-8638 (collection code UJ-13452). One male mature specimen collected at Cherso, Croatia, collected by Marco Curini Galletti 22 September 2013. Photographs of the holotype deposited at the Swedish Museum of Natural History, Stockholm.
Description: Morphologically diverse. Bermuda and Portugal specimens less than 1 mm long, slender, specimens from most other localities comparatively large and bulky; Portuguese specimen was immature. Opaque due to abundance of epidermal glands and oil droplets. Anterior rounded and comparatively slender; mature specimens broaden considerably at level of the statocyst; increase in width more gradual in immature specimens. Posterior tip rounded, possibly with adhesive structure. Frontal glands prominent, opening central, secretion in long bundles; in Bermudian and Portuguese specimens distal base at U35 bulbous, in Adriatic specimens club-shaped (Fig. 8h). Paired testes follicular and lateral; vasa deferentia converge to fsv anterior of mco (Fig. 8g). Male opening in the posterior surrounded by rosette of small glands. Paired ovaries lateral in posterior half or three quarters of body length, oocytes small.
Remarks: This species is supported by the SSU and LSU genes but not by the higher resolving Histone 3 gene. This species is the only truly cosmopolitan species in this dataset. The Histone 3 gene splits the species in geographic clades. This makes the existence of another, not yet resolved species complex, possible. Population genetic studies of this species with increased geographic sampling and sample size might give valuable insights into the dispersal capabilities of interstitial meiofauna.
Etymology: variabilis = latin for variable, as this species is morphologically and genetically changeable throughout its known range. Figure S1 Saturation plot for the H3 gene across the whole dataset. Plotted are the uncorrected p-distances versus the phylogenetic distances between pairs of sequences. The level distribution of the points indicates saturation. (TIFF) Figure S2 Best ML tree calculated with RAxML of the LSU rRNA dataset with bootstrap support plotted on the branches. Putative species with binomial names are formally described in the present study, those with abbreviations represent candidate species. The branch colours correspond to partitions for BP&P analyses, green indicates the Nemertinoides group, red the mainly European Sterreria subgroup and blue the extra-European Sterreria species; orange species have not been validated with BP&P. (TIFF) Figure S3 Majority rule consensus tree estimated with MrBayes of the LSU rRNA dataset with Bayesian posterior probabilities plotted on the nodes. Putative species with binomial names are formally described in the present study, those with abbreviations represent candidate species. The branch colours correspond to partitions for BP&P analyses, green indicates the Nemertinoides group, red the mainly European Sterreria subgroup and blue the extra-European Sterreria species; orange species have not been validated with BP&P. (TIFF) Figure S4 Best ML tree calculated with RAxML of the SSU rRNA dataset with bootstrap support plotted on the nodes. Putative species with binomial names are formally described in the present study, those with abbreviations represent candidate species. The branch colours correspond to partitions for BP&P analyses, green indicates the Nemertinoides group, red the mainly European Sterreria subgroup and blue the extra-European