To refer to this article use this url:

Contributions to Zoology, 85 (2) – 2016

Meiofaunal cryptic species challenge species delimitation: the case of the Monocelis lineata (Platyhelminthes: Proseriata) species complex

Fabio Scarpa1, Piero Cossu1, Tiziana Lai1, Daria Sanna1, Marco Curini-Galletti1, Marco Casu2

1.  Dipartimento di Scienze della Natura e del Territorio – Sezione di Zoologia, Archeozoologia e Genetica, Università di Sassari, Via F. Muroni 25, 07100 Sassari, Italy

2.  Dipartimento di Scienze della Natura e del Territorio – Sezione di Zoologia, Archeozoologia e Genetica, Università di Sassari, Via F. Muroni 25, 07100 Sassari, Italy. E-mail:

Keywords: integrative taxonomic approach,meiofauna conservation,molecular dating,molecular phylogeny,molecular taxonomy,rDNA 18S-28S.


Given the pending biodiversity crisis, species delimitation is a critically important task in conservation biology, but its efficacy based on single lines of evidence has been questioned as it may not accurately reflect species limits and relationships. Hence, the use of multiple lines of evidence has been portrayed as a means to overcome identification issues arising from gene/species tree discordance, morphological convergence or recent adaptive radiations. Here, the integrative taxonomic approach has been used to address the study of the Monocelis lineata species complex. The taxonomic resolution of the complex is challenging, as the species lacks sclerotised copulatory structures, which as a rule of thumb aid identification in Proseriata. Eighteen populations, which encompass most of the geographic range of the complex, were studied using morphology, karyology, crossbreeding experiments and molecular analysis. These different markers provided evidence of four (karyology) to eight (morphology) discrete entities, whereas crossings showed various degrees of intersterility among the tested populations. Molecular species delimitation revealed a different number of candidate species, spanning from five (ABGD and K/θ) to 11 (GMYC). Such incongruences reflect the multifaceted evolutionary history of M. lineata s.l. and hamper the full taxonomic resolution of the complex. However, two candidate species were consistently validated by all of the markers and are described as new species: Monocelis algicola nov. sp. and M. exquisita nov. sp. The latter species appear to have a restricted distribution, and the possibility that meiofaunal taxa may be of conservation concern is discussed.


A recent study of the magnitude of marine biodiversity downscaled the global number of species to the lower end of previous estimates (Appeltans et al., 2012). This result has been challenged on the basis that, in many taxa, the actual number of cryptic species may be underestimated (Mora et al., 2013; Adams et al., 2014). The lack of resolution of cryptic complexes, whose species may have limited geographic ranges and unique ecological requirements (Cox and Moore, 2005), may indeed hinder our chances of determining both the actual level of biodiversity and the extinction risks of the individual components of a cryptic complex.

Molecular techniques are the preferred tool for the resolution of cryptic complexes and, although they have been applied in a minority of instances (see Jörger and Schrödl, 2013 and references therein), nowadays their use is increasing following the implementation of new specific software for automated species delimitation (see Fontaneto et al., 2015). This is particularly true for obscure, minute marine meiofaunal organisms, for which even the description of morphologically distinct species still lags given the tremendous heterogeneity and species-richness of this group. For instance, a recent study of two nominal species of Nemertodermatida, which revealed the existence of 20 supported species (Meyer-Wachsmuth et al., 2014), is of particular interest as it points to the existence of hyper-cryptic complexes in morphologically simple, ‘worm-like’ taxa that lack clear diagnostic features.

Among meiofaunal organisms, the marine microturbellarians (Platyhelminthes) are a group whose contribution to marine biodiversity has yet to be fully assessed (Appeltans et al., 2012). They may be tremendously common and species-rich in littoral marine habitats (Martens and Schockaert, 1986), ranking second in abundance after Nematoda (Fonseca et al., 2010), as evidenced by environmental meta­genomics (Fonseca et al., 2014). This group includes known examples of hyper-cryptic complexes; for instance, in the supposedly cosmopolitan Gyratrix hermaphroditus Ehrenberg, 1831, studies of karyotype and fine morphology have revealed the existence of seven cryptic species in a small portion of a tidal creek at Darwin (NT, Australia) (Curini-Galletti and Puccinelli, 1990). The paucity of diagnostic morphological characters is particularly felt in a taxon of marine microturbellaria, the Monocelididae Monocelidinae (Proseriata), where the copulatory organ of the simplex type (Litvaitis et al., 1996) provides very limited taxonomic information, and species detection and diagnoses increasingly rely on nucleotide-based information (Casu et al., 2009). The case of Monocelis lineata O.F. Müller, 1774 is exemplary as the nominal species has an unusually wide distribution for an interstitial microturbellarian, ranging from both coasts of the North Atlantic through the entire Mediterranean and the Black Sea, and a similarly unusually wide ecological tolerance, occurring in brackish to marine habitats on any type of substrate (Ax, 1956). In stark contrast, most Mediterranean species of Proseriata have ranges limited to single sectors of the basin, and either occur in brackish-water or fully marine habitats (Martens and Curini-Galletti, 1994; Delogu and Curini-Galletti, 2009; Casu et al., 2014). A previous study carried out by means of allozyme electrophoresis (Casu and Curini-Galletti, 2004) on 15 populations morphologically attributable to M. lineata, based on the general morphology of the reproductive structures, suggested the presence of a complex composed of at least six sibling species in Europe alone: five species from the Mediterranean [three (siblings A, B, and E) with a pigmented eye-shield from brackish areas, one unpigmented (sibling C) from sandy marine habitats, and one (sibling D) with a pigmented eye-shield from algae] and one (sibling F) from the Atlantic. However, the study ultimately had no impact on taxonomy due to the lack of more detailed information about markers other than the allozymatic patterns.

This study aimed to provide further insights into the species complex of M. lineata using a multifaceted approach. In fact, since the publication of the paper by Casu and Curini-Galletti (2004), there has been an emerging consensus that an integrative taxonomic approach (Padial et al., 2010; Schlick-Steiner et al., 2010 and references therein), which considers species boundaries from multiple complementary perspectives, is the best way to overcome the potential caveats of any species delimitation method (Dayrat, 2005; Will et al., 2005; Padial et al., 2010).

For the purposes of our study, we planned a workflow that established the six siblings, hereafter treated as OTUs – operational taxonomic units (sensu Sokal and Sneath, 1963; Sneath and Sokal, 1973, identified by Casu and Curini-Galletti (2004) as a taxonomic null hypothesis to be tested against different datasets. The new data sources included information on karyotype, morphology, reproductive biology and genetics (DNA gene sequences). For the genetic approach, we chose to use two nuclear ribosomal markers, the complete nuclear small subunit rDNA (18S) gene and the partial nuclear large subunit rDNA (28S) fragment (spanning the D1-D6 variable domains). These markers have been chosen because of their potential in depicting the genetic variability at inter-specific level, combined to their low level of intra-specific variation, that make them a powerful tool for molecular phylogeny and taxonomy (see Litvaitis et al., 1996; Littlewood et al., 2000; Curini et al., 2010; Casu et al., 2011a; 2014). Furthermore, 18S and 28S genes represent the most complete molecular dataset available for Proseriata in GenBank. Molecular datasets were analysed using automated species delimitation methods, two of which, GMYC (Pons et al., 2006), and PTP/bPTP (Zhang et al., 2013), are based on the PSC (phylogenetic species concept); one, ABGD (Puillandre et al., 2012) is based on genetic distances, and one, K/θ method (Birky et al., 2010) is based on the 4× rule criterion. Finally, we estimated the divergence time among taxa to obtain a better understanding of the evolutionary pathways (see Lemey and Posada, 2009).

Material and methods

Taxon sampling

We utilised specimens from populations of the Mono­celis lineata complex evinced as distinct OTUs by Casu and Curini-Galletti (2004) (and unpubl. data) (see Appendix 1 for details). To gain a wider taxonomic landscape, we designed a molecular dataset with the available species of the genus Monocelis: M. longiceps (Dugés, 1830), M. longistyla Martens & Curini-Galletti, 1987, M. fusca Örsted, 1843 and other species belonging to the subfamily Monocelidinae: Minona ileanae Curini-Galletti, 1997; Pseudomonocelis ophio­cephala (Schmidt, 1861); P. occidentalis Curini-Galletti, Casu & Lai, 2011; P. orientalis Curini-Galletti, Casu & Lai, 2011; P. agilis (Schultze, 1851); P. cetinae Meixner, 1943; P. cf. cavernicola Schockaert & Martens, 1987; P. paupercula Curini-Galletti, Casu & Lai, 2011. Archiloa rivularis de Beauchamp, 1910, a Mono­celididae Duplomonocelidinae (sensu Litvaitis et al., 1996), and Archimonocelis staresoi Martens & Curini-Galletti, 1993 (Proseriata: Archimonocelidae) were chosen as outgroups.

Samplings and isolation of specimens

Samples were collected manually by scooping up the superficial layer of sediment, or by scratching algae from the upper infralittoral. Permits for samplings in the National Park of La Maddalena (Italy) were obtained from the Director of the Park, Dr. M. Gargiulo. No specific permits were required for other sites, which were not privately owned or protected. Given the size of the organisms, the samples collected were small (in most instances, less than 2 l of sediments). Extraction of the animals from the sediment and from algae was accomplished using MgCl2 decantation (see Casu et al., 2011a and reference therein). Each individual was first studied alive by slight squeezing under the cover slip. Specimens were then retrieved and processed for further analyses.

Morphological and karyological analysis

Twenty-two to 33 fully mature specimens per population, after relaxation in an isotonic MgCl2 solution, were fixed in Bouin’s fluid and embedded in Paraplast at 56°. Serial sections were cut at 4 μm thickness, stained in Mayer’s haematoxylin and eosin, and mounted in Depex.

Biometric analyses took into account the following morphological parameters: (1) total body length; (2) length and (3) width of the copulatory organ; (4) length of the penis papilla; (5) thickness of the proximal muscle sheath of the copulatory bulb; (6,7) thickness of the (6) muscular and (7) prostatic components of the penis papilla; (8) distance between vaginal pore and male pore; (9) distance between male pore and female pore. All measurements were performed by the same person (MCG) on the same microscope (Leitz DM-RB).

In order to obtain karyometrical information, 15 to 20 specimens per population were placed in a 2% colchicine solution for 6 h. Specimens were then put on a slide, and kept for 5 min in acetic acid 5%. After removal of the solution, specimens were stained with lactic acetic orcein for about 7 min, covered with a coverslip, and strongly squeezed. This staining technique, which can be carried out in field conditions, as often happened in the course of the research, and with minimal equipment, does not allow removal of tissues. For this reason, chromosome plates rarely lie perfectly at the same focal plane, and, in many instances, need to be studied through slight movements of the micrometric focusing knob, and were unsuitable for photomicrography. Measurements were thus obtained from metaphase plates drawn utilising a camera lucida. The measured parameters were: relative length (r.l.= length of chromosome × 100/total length of haploid genome) and centromeric index (c.i.= length of short arm × 100/length of entire chromosome).

To test whether the karyological and morphological characters may yield taxonomic information, unifactorial analysis of variance (ANOVA) was performed on each variable. Cochran’s C-test was used to test the assumption of homogeneity of variance. Following significant effects in the ANOVA (P < 0.01), a posteriori analyses were performed [Student–Newman–Keuls (SNK) test], in order to detect homogeneous groupings of populations at a P = 0.05 (Underwood, 1997).

In order to assign individuals to clusters we followed the approach described in Edwards and Knowles (2014), which combines multivariate and clustering techniques, with modifications. The method allows to combine data of different origin (e.g., morphological, genetic, ecological, etc.) using multivariate non-metric multidimensional scaling (nMDS) to obtain a new set of reduced variables (the nMDS dimensions) from each independent dataset. Such standardized datasets are then analysed separately or jointly using Gaussian clustering; in the second case reduced nMDS variables from each dataset are merged in a unique matrix. However, in this study it was not possible to combine standardized dimensional datasets in a single matrix, since morphological and karyological as well as genetic analyses could not be performed on the same set of individuals, due to the small size of the flatworms. All analyses and packages reported above are implemented in the R 2.15.3 statistical environment (R Core Team, 2013).

Crossbreeding experiments

Cultures of specimens were maintained in the laboratory at 18 ± 0.5 °C and 37‰ salinity and fed weekly with crushed Sphaeroma sp. (Isopoda). Reproductive biology experiments were carried out on one population for each of the OTUs proposed by Casu and Curini-Galletti (2004) (and unpubl. data): ALo (OTU A), CSo (OTU B), PPx (OTU C), CRo (OTU D), PUo (OTU E), KFo (OTU F). They all showed good survival and remarkable longevity in the laboratory. Two sets of experiments were carried out:

a) Intra-population crossings:

Experiments were carried out by isolating pairs of newborns in 20 ml containers filled with seawater, under the culture conditions described above. The offspring produced weekly by pairs were counted and relocated. Fifteen replicates (i.e., pairs) were set for each of the six populations. To exclude bias due to ageing of specimens, which results in sparse, episodic events of reproduction (unpubl. data), the experiments were carried out for a 20 weeks period. In case one or both members of the pair died before 20 weeks from the onset of the experiment, the replicate was considered as lost, and a new pair was set as replacement.

b) Inter-population crossings:

Experiment was set as above, but the pair of juveniles belonged to different populations. All possible combinations were tested (i.e., ALo × CSo; ALo × PPx; ALo × CRo; ALo × PUo; ALo × KFo; CSo × PPx; CSo × CRo; CSo × PUo; CSo × KFo; PPx × CRo; PPx × PUo; PPx × KFo; CRo × PUo; CRo × KFo; PUo × KFo), with 15 replicates per combination (for details about abbreviation codes see Appendix 1).

ANOVAs were used to analyse differences in reproductive outputs (i.e., the number of offspring produced) among pairs of specimens of different populations in the 20 weeks period. Cochran’s C-test was used to test the assumption of homogeneity of variance. The SNK test was used a posteriori to compare means.

DNA extraction, amplification and sequencing

Genomic DNA was extracted using the Macherey-Nagel NucleoSpin Tissue (MACHEREY-NAGEL GmbH & Co. KG) according to the supplier’s instructions. After extraction, DNA was stored as a solution at 4 °C. Complete 18S and partial 28S (spanning variable domains D1-D6) were analysed for a total of 56 specimens, 45 of which newly sequenced by us, and 11 taken from GenBank (for details about specimens and sampling localities see Appendix 1).

PCRs for 18S and 28S D1-D6 regions were carried out using the following primers: 18S: A (forward) GCG AAT GGC TCA TTA AAT CAG, and B (reverse) CTT GTT ACG ACT TTT ACT TCC (Littlewood and Olson, 2001); 28S D1-D6: LSU5 (forward) TAG GTC GAC CCG CTG AAY TTA AGC A, and LSUD6-3 (reverse) GGA ACC CTT CTC CAC TTC AGT C (Littlewood et al., 2000). PCRs were carried out in a total volume of 25 μl containing 5 ng/μl of total genomic DNA on average, 1.0 U of Taq DNA Polymerase (Euroclone), 1× reaction buffer, 3.5 mM of MgCl2 , 0.32 μM of each primer, and 200 μM of each dNTP. PCR amplifications were performed in a MJ PTC 200 Thermal Cycler (Biorad) programmed as follows: 1 cycle of 2 min at 94° C, 35 cycles of 1 min at 94° C, 1 min at 54° C (18S / 28S D1-D6 primers’ annealing temperature), and 1 min and 30 s at 72° C. At the end, a post-treatment for 5 min at 72° C and a final cooling at 4° C were carried out. Both positive and negative controls were used to test the effectiveness of the PCR reagents, and the absence of possible contaminations.

Electrophoresis was carried out on 2% agarose gels, prepared using 1× SBA buffer (sodium boric acid, pH 8.2) and stained with a 1 μl/20 ml ethidium bromide solution. PCR products were purified by ExoSAP-IT (USB Corporation) and sequenced for both forward and reverse 18S and 28S D1-D6 strands (by means of the same primers used for PCR), using an external sequencing core service (Macrogen Europe).

Phylogenetic analysis and dating

The 18S and 28S D1-D6 sequences were aligned separately using the algorithm Q-INS-I implemented in Mafft 6.903 (Katoh and Toh, 2008), which is appropriate for non-coding RNA as it considers RNA secondary structure. The best probabilistic model of sequence evolution was determined independently for each gene using jModeltest 2.1.1 (Posada, 2008), with a maximum likelihood optimised search, and both the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC). The model GTR+G (Tavaré, 1986) both for 18S and 28S D1-D6 was selected by either criterion as the best-fitting model to both sequences.

Phylogenetic relationships among taxa were investigated using both Maximum Likelihood (ML) and Bayesian Inference (BI) on the combined 18S and 28S D1-D6 sequences.

ML analyses were conducted using the software RAxMLGUI version 1.3 (Silvestro and Michalak, 2012) setting the analysis option to ‘ML+ thorough bootstrap’, which consists in a thorough bootstrap analysis followed by a maximum likelihood search. Then, the bootstrap support values are drawn on the most likely tree. For each data partition, separate GTR+G models were used. Analysis was carried out with 10 runs and 1000 bootstrapping replicates. The consensus tree was visualized by means of the software FigTree 1.4.0 (

BI was carried out using the software package Beast 1.7.4 (Drummond and Rambaut, 2007). This software was also used to obtain both the ultrametric tree and the chronogram for estimating the divergence time from the most recent common ancestor (TMRCA) for each clade. Site parameters (Substitution Model = GTR; Bases Frequencies = Estimated; Site Heterogeneity Model = Gamma; Number of Gamma Categories = 4) have been set according to the model selected by jModeltest. The lognormal uncorrelated relaxed molecular clock model was chosen because it assumes independent mutation rates on different branches. For the tree prior, the Yule prior process to the speciation model was applied.

Priors for model parameters and statistics have been set using a normal distribution, with lower and upper values set according to Scarpa et al. (2015). Operator parameters have been set following the instructions on the user manual. In order to obtain an Effective Sample Size (ESS) greater than 200 for all of the parameters, a run of 200,000,000 generations was performed, sampling a tree every 20,000 generations. Run was executed by means of the Cipres Phylogenetic Portal (Miller et al., 2010). Resulting log files were examined using the software Tracer 1.6 (Rambaut and Drummond, 2009). TreeAnnotator (part of the BEAST package) was used for drawing the chronogram, applying a Burn-in period of 10%. The chronogram has been visualised by means of the software FigTree (

Alignments and Bayesian tree-files are deposited and available in Treebase with the accession number TB2: S17581.

Species delimitation

We used four different methods of species delimitation. First we applied two methods inspired to the PSC (Phylogenetic Species Concept), which use the species tree to delimit species: the GMYC (Generalized Mixed Yule Coalescent) method by Pons et al. (2006), and the PTP (Poisson Tree Processes) model and its Bayesian implementation, the bPTP (Zhang et al., 2013).

The GMYC model tests for a significant shift in the branching rate along an ultrametric species tree. The analysis was performed on the ultrametric species tree resulting from the Bayesian analyses Fig. 1() by means of the SPLITS (SPecies LImits by Threshold Statistics) package (Ezard et al., 2009) implemented in the R statistical environment (available at As recommended by the user manual, species entities were identified by means of the single threshold option, which uses a single threshold to specify the transition from between- to within-species branching.


Fig. 1. Chronogram obtained by Bayesian Inference (BI) showing interrelationships of the species based on combined 18S + 28S D1-D6. The branch length is proportional to the number of substitutions per site. Node supports are indicated by posterior probabilities for BI, and bootstraps for Maximum likelihood (ML). The grid indicates the mean divergence time (in myr). A synthetic stratigraphic chart, in which geological epochs have been placed according to Cohen et al. (2013), is reported below the time bar.

The PTP and bPTP models use the number of substitutions to assess the speciation rate. Species delimitation was performed by means of the bPTP web server (available at on the Bayesian phylogenetic species tree (Fig. 1(), using default options and 500,000 MCMC generations. In order to test the reliability of results, each run was checked for convergence by visualising the likelihood plot: if convergence occurred, the chain should stay at high likelihood locations most of the time during the run.

We also applied the ABGD (Automatic Barcode Gap Discovery) (Puillandre et al., 2012) method on the combined dataset (18S + 28S D1-D6) using K2P genetic distances (Kimura, 1980). This method does not consider phylogenetic relationships within the dataset, and works on sequences, detecting the barcode gap as the first significant gap beyond this limit, and using it to partition the data (Puillandre et al., 2012). Species were assessed by means of the ABGD online tool (available at using the default settings. The correct species estimate was selected, as suggested by Puillandre et al. (2012), using the gene specific priors for maximum divergence of intraspecific diversity, corresponding to P = 0.001.

The last method applied is the K/θ (Birky et al., 2010). This method is based on a speciation model relying upon the so-called 4x rule criterion. According to this criterion, if the ratio between the mean pairwise sequence difference between a pair of clades (K), and the mean pairwise sequence difference within a clade (θ) is greater than 4, clades are samples from different species (Birky, 2013). Thus, the K/θ method estimates the probability that two samples are from different evolutionary species (Birky et al., 2010). For this test the topology obtained in present study (see Fig. 1() has been used. K/θ has been applied on a subset composed by specimens belonging to the M. lineata complex sensu stricto, OTUs D and E, and M. fusca. On this subset the uncorrected (p-distance) pairwise genetic distances have been calculated by means of the software Mega 6.06 (Tamura et al., 2011) with 1,000 bootstrap replications. For estimating the ratio K/θ, formulas corrected for sexual organism have been used as exposed in Birky (2013). Clades that showed values of K/θ ≥ 4 should be considered as entities with a 95% probability of having an independent evolutionary history. In order to be as conservative as possible, for each ratio, the upper values of θ have been used.


Morphological data

Values of the nine morphological parameters are shown in S1. ANOVAs revealed significant differences in the sample for most of the variables tested, but groupings, as detected by SNK test (S2), were not consistent among variables.

Stress value of the nMDS ordination was 8.87%, and the plot of the first two nMDS dimensions did not evidence a clear separation between individuals belonging to different OTUs, with the exception of CRo (OTU D) (Fig. 2(). In addition, the eight clusters identified by the Gaussian clustering (S3A) of the standardized morphological nMDS dimensions showed a poor congruence either with the six OTUs based upon the allozyme data (Casu and Curini-Galletti, 2004), or with geography, as well as presence/absence of pigmented shield. Indeed, with respect to the six OTUs the classification error was as high as 66%. Nonetheless, two Gaussian clusters matched tightly OTU D (CRo) and OTU E (PUo), respectively (S3A). Using a probability threshold of 0.95 to assign individuals to a species cluster, only two specimens from PUo could not be assigned to any cluster. In the remaining populations the percentage of unassigned individuals ranged from 12% (PPx) to 76% (HEo), with the exception of CSo where all specimens were assigned to the same cluster with a probability ≥ 0.95. Interestingly, successfully assigned Mediterranean and Atlantic individuals were never assigned to the same cluster.


Fig. 2. Non-Metric Multidimensional Scaling (nMDS) ordination of morphological data. The plot shows the relationships among individuals along the first two dimensions. Mediterranean and Atlantic locations are indicated with circles and triangles, respectively.


Fig. 3. Non-Metric Multidimensional Scaling (nMDS) ordination of karyological data. The plot shows the relationships among individuals along the first two dimensions. Mediterranean and Atlantic locations are indicated with circles and triangles, respectively.

Karyological data

Values of the six karyological parameters are shown in S4. ANOVAs performed on karyological variables revealed significant differences among populations, but SNK failed to detect discrete groupings of populations (S5), apart from the centromeric index of Chromosome III, where Mediterranean and Atlantic populations (with the exception of CRo) were assigned to different groupings. Furthermore, CRo was uniquely characterised for centromeric index of Chromosome II, and PUo for relative lengths of Chromosomes I and III, and centromeric index of Chromosome III.

The plot of the first two nMDS dimensions separated CRo and PUo (OTUs D and E, respectively) better than morphological data (Figs 2(, 3(); the stress value of the ordination was slightly lower as well (7.24%). Gaussian clustering of standardized dimensional karyological data detected four clusters (S3B). The classification error rate with respect to the OTUs was less than that corresponding to morphological data (30%). Two Gaussian clusters corresponded almost perfectly to CRo and PUo, respectively (S3B), and nearly all individuals were assigned to the expected cluster with a probability of ≥ 95%. Once again two specimens from PUo could not be assigned to any cluster as their probability did not reach the assumed threshold; conversely an individual from COo was assigned to the same cluster as individuals from PUo. The remaining individuals that could be assigned to a cluster with a probability ≥ 95% were subdivided into two clusters that roughly resembled their Mediterranean or Atlantic origin. Only one Atlantic (from FEo) and eight Mediterranean specimens (four from PPx and four from CHx) were misassigned. Furthermore, with the exceptions of CRo and TJx (where all individuals were assigned), the proportion of unclassified individuals ranged from about 7% (ARx) up to 80% (CSo and PLo).

Crossbreeding experiments

Intra-populational crossing experiments revealed significant differences in the number of offsprings produced by pairs during the 20 weeks period of the test (S6) (Cochran’s C Test= 0.2792, P>0.05; F=57.21, P<<0.001); SNK test: PPx=CRo<<ALo=KFo=PUo<< CSo).

On the contrary, no offsprings were produced in any of the inter-populational crossings performed, with the unique exception of the test involving pairs derived from specimens from CSo and ALo, which produced 1.7437 ± 0.8681 offsprings per week.

Phylogeny and dating

After the alignment, 1656 and 1596 bp-long sequences were obtained for the 18S and 28S D1-D6 regions, respectively (see Appendix 1 for the GenBank accession numbers). Because both the ML and BI results converged on the same topology, we reported only the Bayesian tree (Fig. 1(). All of the main nodes in the tree are highly supported by both bootstrap values and posterior probability.

Populations attributed to the Monocelis lineata species complex were not monophyletic as Monocelis fusca (MFU) is nested among them. In particular, MFU is the sister taxon of PUo, whereas CRo shows a sister-taxon relationship with the clade including the rest of the Mediterranean and Atlantic populations. The remaining Mediterranean samples constitute a monophyletic group split into two reciprocally monophyletic clades, one composed of individuals from marine habitats and one of individuals from brackish habitats (Appendix 1). Samples from the brackish habitats were further split into two clusters representative of the two brackish taxonomic units retrieved by allozyme electrophoresis analyses (OTUs A and B). Conversely, individuals belonging to the Atlantic populations do not form a monophyletic clade but are paraphyletic because FEo clustered externally to the Mediterranean populations. Remarkably, the tree shows that the Monocelis genus itself is paraphyletic; indeed, Monocelis longiceps (MLS) results as an outgroup of the clade composed by all of the specimens belonging to the genus Pseudomonocelis + Minona ileanae (MIL), although this clustering is very low supported.

The estimated ucld.stdev parameter amounts to 0.977 and 0.691 for 18S and 28S, respectively, indicating that our dataset is clock-like (i.e., it does not show substantial rate heterogeneity among lineages). The divergence times for the whole dataset are reported in Fig. 1(, S7.

Molecular species delimitation

The results for the whole dataset are as follows: i) the GMYC model identified a total of 24 entities (CI = 8-27), 12 of which are represented by singletons and 12 by clusters (CI = 2-13) (likelihood ratio test of the null against the mixed model: 5.4e-05, P < 0.001); ii) the PTP/bPTP method identified a total of 20 entities (CI = 16-25), and iii) the ABGD method, checked at the prior maximal distance (P = 0.001), identified 15 entities.

Regarding the M. lineata complex i) the GMYC model found five entities among the Mediterranean samples and six in the Atlantic; ii) the PTP/bPTP model found a total of seven entities, three from the Mediterranean and four from the Atlantic; iii) the ABGD method found four entities in the Mediterranean and one in the Atlantic, and iv) the K/θ method delimited three and two entities from the Mediterranean and Atlantic, respectively. For details on the species delimitation results, see Fig. 4( and S8, S9, S10, S11.


Fig. 4. Summary of the results obtained by all species delimitation methods applied. The first row refers to the six OTUs obtained in a previous paper (Casu and Curini-Galletti, 2004), here used as null hypothesis. For morphological and karyological datasets, colours and patterns inside the squares correspond to the clusters identified using Gaussian clustering to which specimens have been assigned, respectively. Only those individuals that have been successfully assigned to a cluster above the 95% probability threshold have been reported in the figure. Squares containing different colours or patterns correspond to populations in which individuals were successfully assigned to more than one cluster. Amount of colours or patterns is proportional to the number of individuals assigned to each cluster.


The study of free-living Platyhelminthes can be challenging due to the limited number of characters with taxonomic and phylogenetic weight. However, in some interstitial microturbellarian taxa, such as the Rhabdocoela, the complex morphology of the sclerotised pieces of the copulatory organ (and, in some instances, of the female reproductive system) allows for species discrimination with a certain degree of confidence (see e.g., Artois, 2008). However, these characters are not always present; among the order Proseriata, the Monocelididae Monocelidinae usually lack sclerotised parts in their copulatory organ (see e.g., Litvaitis et al., 1996), substantially reducing the number of characters available for study. This limitation calls for the utilisation of markers other than morphology, particularly molecular markers, which have significantly contributed to the detection of species boundaries (e.g., Casu et al., 2011a; Casu et al., 2014; Leasi and Norenburg, 2014 and references therein).

In our study, we considered the six OTUs identified by Casu and Curini-Galletti (2004) as candidate species (null hypothesis). Multiple lines of evidence (i.e., morphometric and karyometric data; reproductive biology; genetic data from two nuclear loci, analysed using a phylogenetic approach and four different species delimitation methods) were tested on 18 populations, which span most of the range of the species complex. The patterns revealed by these markers were not entirely congruent:

Morphology: Of the eight clusters identified, two matched OTUs D and E (CRo and PUo, respectively). The rest of the Mediterranean and all of the Atlantic specimens were grouped into six clusters with no correspondence to the OTUs (Figs 4(, S3A).

Karyology: Two of the four clusters detected by Gaussian clustering matched OTUs D and E, with the exceptions mentioned above. However, specimens from the populations attributed to the four OTUs (A - C, F) were assigned to only two clusters, which were predominantly but not exclusively based on geography (Mediterranean vs Atlantic) (Figs 4(, S3B).

Reproductive biology: The experiment indicated a barrier to crossbreeding among the studied populations, with the exception of specimens of OTUs A (ALo) and B (CSo), which showed reduced fertility compared to crossings involving specimens from the parental populations. On the other hand, PLo (also belonging to OTU B), did not crossbreed with other OTUs (S6). Distinct allopatric lineages may show reproductive compatibility (de Queiroz, 2005), and the results do not necessarily disprove the null hypothesis, which appears to be supported by the presence of reproductive barriers among the rest of the study populations.

Molecular analysis: Our phylogenetic analysis indicated that Monocelis fusca (MFU), which was initially assumed to be the sister taxon to our study case, clustered within the M. lineata complex. The status of M. fusca as a valid species has never been questioned: among other features, it is characterised by the presence of a copulatory stylet variable in length (see Graff, 1913). The other Monocelis species with a copulatory stylet in our database, M. longistyla (MLA) appears to be only remotely related to M. fusca (Fig. 1(), which raises questions about how a stylet may be acquired/lost in Monocelididae. However, this is beyond the scope of our paper.

Within the M. lineata complex (henceforth excluding M. fusca), the results of the molecular-based species delimitation methods were not consistent as they detected different numbers of entities (Fig. 4(). For instance, within the Mediterranean populations, GMYC merged the brackish OTUs A and B and split the marine populations (OTU C) into two entities, from the western (PPx) and eastern Mediterranean (CHx). In contrast, the K/θ indicated no separation among OTUs A, B and C. Even more heterogeneous are the results obtained for the Atlantic populations (Fig. 4() from which we found one (ABGD) to six (GMYC) different entities. Although we cannot rule out that the different numbers of entities identified by the four methods may reflect the different modelling assumptions (Carstens et al., 2013), the recent divergence time may account for the discrepancy in the number of entities obtained by the species delimitation methods (Fig. 4(), as detecting species boundaries tends to be less accurate when divergence times are low (Puillandre et al., 2012; Birky, 2013; Fujisawa and Barraclough, 2013). For instance, our results dated the basal node of the A, B and C OTUs to approximately 800,000 years ago (Fig. 1(, S7) and set the divergence time between PPx and CHx at only 530,000 years ago (Fig. 1(, S7). In addition, within the Atlantic entities, divergences dated from approximately 700,000 to 330,000 years ago (Fig. 1(, S7).

In this context, it is surprising that GMYC, one of the most widely used methods (see Tang et al., 2012 and references therein), identified quite a few Atlantic populations as distinct entities, i.e., FEo, HEo, RSx, TJx, and KQo (Fig. 4(). These populations are located in geographic areas recognised as Quaternary age refugia for M. lineata (Casu et al., 2011b). Thus, the complex phylogeography of Monocelis lineata s.s. in the Atlantic might have clouded the pattern, as intraspecific divergence can accumulate over evolutionary time in geographically isolated populations (Norris and Hull, 2012). It should be noted that the rate of genome evolution may vary significantly according to species characteristics, ecological factors, aspects of evolutionary history and life history traits (Thomas et al., 2010 and references therein). Accordingly, poor dispersal capabilities, as shown by the Proseriata (Curini-Galletti et al., 2012) may promote adaptation to local conditions and enhance intraspecific genetic and phenotypic variation.

Notwithstanding the largely unresolved taxonomic scenario, two populations, CRo and PUo (OTUs D and E, respectively), were consistently validated by multiple lines of evidence (see Appendix 2; see also Tables 1, 2 for molecular pure diagnostic characters) and represent candidate species belonging to the CCS category (confirmed candidate species), i.e., entities that can be delimited by molecular data and are also supported by other data but have not yet been formally described and named (see Vieites et al., 2009). There may be few objections to attributing the rank of species to these highly characterised lineages. Indeed, in addition to the molecular data, the two populations are uniquely characterised by morphological and/or karyological data (see Appendix 2). This de facto outcome reduces the M. lineata complex to OTUs A, B, C, and F of the primary hypothesis, with the two new species and M. fusca as sister taxa, and stresses the importance of extensive taxonomic and geographic sampling.


Table 1. Molecular pure diagnostic characters based on the 18S and 28S genes for CRo. Nucleotide positions refer to CRo sequences (GenBank accession number: KR364646-50 for the 18S gene and KR364691-95 for the 28S D1-D6 fragment).


Table 2. Molecular pure diagnostic characters based on the 18S and 28S genes for PUo. Nucleotide positions refer to PUo sequences (GenBank accession number: KR364641-45 for the 18S gene and KR364686-90 for the 28S D1-D6 fragment).

Naming the species

Although integrative taxonomy gives priority to species delineation and not to the creation of new names (Dayrat, 2005), the importance of naming species that result from the validation of the primary species hypotheses should not be underestimated. Without a formal description, in fact, newly discovered species may not be documented and adequately vouchered, and reproducibility can thus be hindered and confusion can result from the different numbering systems of the OTUs (Jörgen and Schrödl, 2013).

The synonymy of Monocelis lineata is vast (see Graff, 1913), and the first nomenclatural problem is determining to which of the populations the nominal taxon applies. The taxon was introduced by O.F. Müller (1774: pp 60, 61) for specimens found ‘in littore maris Balthici’ with a pigmented eye-shield. At present, M. lineata is known to occur throughout most of the Baltic, from the Øresund to Finland (Luther, 1960; Karling, 1974). Baltic specimens have a vertically elongate, elliptical copulatory organ and a vagina-male pore distance that is markedly longer than the male pore-female pore distance (Luther, 1960). In our sample, the pigmented population that is closer to the type locality is HEo from the Øresund. Indeed, this population shows a muscular, ovoid-elongate copulatory organ; the pore ratio is 3.68: 2.12: 1 (based on 15 specimens; see Appendix 2 for the explanation of the ratio) (Fig. 5I-M(). HEo specimens have a karyotype with isobrachial chromosomes that are almost even in size (Fig. 5M(, S4). To the best of our knowledge, this population may be taken as representative of the nominal species. Most of the taxa synonymised with M. lineata are based on populations from the NE Atlantic Ocean (Graff, 1913), and the original, extremely poor descriptions and the lack of museum specimens for comparison make their classification problematic if not impossible (Karling, 1966) and should best be considered to be nomina dubia.

In any case, the taxonomic classification proposed for the northern European populations does not appear applicable to the two Mediterranean populations. They are distinct from any of the known Atlantic populations for almost any marker considered and are characterised by univocal morphological, karyological, and molecular markers (this paper). We therefore describe PUo (OTU E) and CRo (OTU D) as new species, namely Monocelis exquisita nov. sp. and Monocelis algicola nov. sp., respectively (see Appendix 2).

A conservation issue?

A shared trait of the new species is their apparent restricted distribution; both are known from only one station. Monocelis algicola nov. sp. is the only Mediterranean Monocelis species occurring on algae. Since the habitat has not been thoroughly sampled elsewhere for Proseriata, the species may be more widespread than presently known.

Monocelis exquisita nov. sp. may prove to be a different case. It is only known from a shallow brackish water habitat in Northern Sardinia (Porto Puddu). Specimens were first sampled in 2002 (Casu and Curini-Galletti, 2004), and at that time, the species was present throughout the inlet and in an adjacent creek mouth. Afterwards, the area underwent massive modifications due to the construction of a pier and a marina (S12) and a real estate development. Extensive sampling performed in June 2014 revealed that the species had possibly disappeared from most of the stations where it was common in 2002, and its presence was only confirmed in the unaffected creek (S12). Closely related species of the genus Pseudomonocelis Meixner, 1943 have strict sediment preference, and populations have disappeared after natural or human mediated sedimentary imbalance (Cognetti and Curini-Galletti, 1993; Casu and Curini-Galletti, 2006). The apparent disappearance of the inlet population of M. exquisita nov. sp. may thus be linked to the observed silting of the sediments of the area, presumably caused by reduced exchange with the sea following the construction of the pier.

Brackish water habitats, and especially brackish microhabitats, are numerous throughout Sardinia, and it would be pretentious to assume that we managed to sample all of the habitats suitable for the species. However, over the years, extensive sampling along most of the Sardinian coast (our own unpubl. data) failed to detect the species. Instead, we consistently found one or the other of OTUs A and B. Brackish areas were not as thoroughly investigated for Proseriata outside of Sardinia, but we have conducted numerous sampling campaigns in many parts of the northern and central Mediterranean (mostly unpubl. data). In most of the brackish areas of the Mediterranean, one or the other of OTUs A and B were abundantly found, but in none of them were specimens of M. exquisita present. The demise of the species in its only known station is thus particularly worrying.

The species meets the IUCN criteria of Critically Endangered (CR), as the species is: “... known to exist at only a single location.” (point V, B2.a), and the area, extant and quality of habitat is in “continuing decline, observed, inferred or projected ...” (point V, B2.b iii) (IUCN Red List Categories and Criteria. Version 3.1, 2001).


Discrepancies among markers and methods may contribute to the elucidation of their limits. Indeed, although genetic data have recently became the primary workhorse used for species delimitation, an uncritical adoption of any species delimitation method in this case would have resulted in a taxonomic resolution that might have not reflected the evolutionary history and speciation processes of the M. lineata complex.

Purely genetic methods can inadequately describe diversity and require significant amounts of genetic data to delimit rapid radiation (Edwards and Knowles, 2014 and references therein), and evolutionary processes may impact the detection of species using only a single datatype (either genetic or morphological). The integration of information from as many markers as possible, i.e., integrative taxonomy, remains thus the best option for overcoming the limits of a single source of data.

Although we partially resolved the M. lineata complex by describing two new species, there is evidence pointing to the existence of an Atlantic-Mediterranean hyper-cryptic complex, which still needs to be explored. More researches involving more populations across biogeographic breaks, more inclusive cross-breeding experiments, and the adoption of molecular markers better tailored to tackle the rapid radiations that occur along morphological/ecological axes (see Wagner et al., 2013 and references therein) are still needed.

Our research has involved a number of populations and extensive coverage of a geographic area, which is unusual in studies of meiofauna, and a potentially interesting outcome is that these minute organisms, which are usually completely neglected in overviews of marine biodiversity conservation (e.g., Bianchi, 2007), may have restricted distributions (see also Casu and Curini-Galletti, 2006; Casu et al., 2009, 2011a, 2012, 2014; Leasi and Norenburg, 2014), and can be particularly affected by the alteration of habitats, facing the threats of silent extinction.


We acknowledge for partial funding support ASSEMBLE grant agreement no. 227799 for specimens sampled in Roscoff. Valentina Delogu skillfully sectioned the specimens used in this study. This research is part of the PhD thesis of Fabio Scarpa.

Received: 23 April 2015

Revised and accepted: 22 September 2015

Published online: 21 March 2016

Editor: R. Sluys


Adams M, Raadik TA, Burridge CP, Georges A. 2014. Global biodiversity assessment and hyper-cryptic species complexes: more than one species of elephant in the room? Systematic Biology 63: 518-33.

Appeltans W, Ahyong ST, Anderson G, Angel MV, Artois T, Bailly N, Bamber R, Barber A, Bartsch I, Berta A, Błazewicz-Paszkowycz M, Bock P, Boxshall G, Boyko CB, Brandão SN, Bray RA, Bruce NL, Cairns SD, Chan T-Y, Cheng L, Collins AG, Cribb T, Curini-Galletti M, Dahdouh-Guebas F, Davie PJF, Dawson MN, De Clerck O, Decock W, De Grave S, de Voogd NJ, Domning DP, Emig CC, Erséus C, Eschmeyer W, Fauchald K, Fautin DG, Feist SW, Fransen CHJM, Furuya H, Garcia-Alvarez O, Gerken S, Gibson D, Gittenberger A, Gofas S, Gómez-Daglio L, Gordon DP, Guiry MD, Hernandez F, Hoeksema BW, Hopcroft RR, Jaume D, Kirk P, Koedam N, Koenemann S, Kolb JB, Kristensen RM, Kroh A, Lambert G, Lazarus DB, Lemaitre R, Longshaw M, Lowry J, Macpherson E, Madin LP, Mah C, Mapstone G, McLaughlin PA, Mees J, Meland K, Messing CG, Mills CE, Molodtsova TN, Mooi R, Neuhaus B, Ng PKL, Nielsen C, Norenburg J, Opresko DM, Osawa M, Paulay G, Perrin W, Pilger JF, Poore GCB, Pugh P, Read GB, Reimer JD, Rius M, Rocha RM, Saiz-Salinas JI, Scarabino V, Schierwater B, Schmidt-Rhaesa A, Schnabel KE, Schotte M, Schuchert P, Schwabe E, Segers H, Self-Sullivan C, Shenkar N, Siegel V, Sterrer W, Stӧhr S, Swalla B, Tasker ML, Thuesen EV, Timm T, Todaro MA, Turon X, Tyler S, Uetz P, van der Land J, Vanhoorne B, van Ofwegen LP, van Soest RWM, Vanaverbeke J, Walker-Smith G, Walter TC, Warren A, Williams GC, Wilson SP, Costello MJ. 2012. The magnitude of global marine species diversity. Current Biology 22: 2189-2202.

Artois TJ. 2008. Revision of Rogneda Uljanin, 1870 (Rhabdito­phora, Eukalyptorhynchia, Polycystididae) with the description of seven new species. Zoological Journal of the Linnean Society 153: 1-28.

Ax P. 1956. Monographie der Otoplanidae (Turbellaria). Morphologie und Systematik. Abhandlungen der Matematisch-naturwissenschaftlichen Klasse 1. Akademie der Wissenschaften und der Literatur Mainz 13: 499-796.

Bianchi CN. 2007. Biodiversity issues for the forthcoming tropical Mediterranean Sea. Hydrobiologia 580: 7-21.

Birky CW Jr, Adams J, Gemmel M, Perry J. 2010. Using population genetic theory and DNA sequences for species detection and identification in asexual organisms. PLoS One 5: e10609.

Birky CW Jr. 2013. Species detection and identification in sexual organisms using population genetic theory and DNA sequences. PLoS One 8: e52544.

Carstens BC, Pelletier TA, Reid NM, Satler JD. 2013. How to fail at species delimitation. Molecular Ecology 22: 4369-4383.

Casu M, Curini-Galletti M. 2004. Sibling species in interstitial flatworms: a case study using Monocelis lineata (Proseriata: Monocelididae). Marine Biology 145: 669-679.

Casu M, Curini-Galletti M. 2006. Genetic evidence for the existence of cryptic species in the mesopsammic flatworm Pseudomonocelis ophiocephala (Rhabditophora: Proseriata). Biological Journal of the Linnean Society 87: 553-576.

Casu M, Lai T, Sanna D, Cossu P, Curini-Galletti M. 2009. An integrative approach to the taxonomy of the pigmented European Pseudomonocelis Meixner, 1943 (Platyhelminthes: Proseriata). Biological Journal of the Linnean Society 98: 907-922.

Casu M, Cossu P, Sanna D, Lai T, Scarpa F, Curini-Galletti M. 2011a. A reappraisal of the monophyly of the genus Pseudo­monocelis Meixner, 1943 (Platyhelminthes: Proseriata), with the description of a new species from the Mediterranean. Zootaxa 3011: 59-68.

Casu M, Sanna D, Cossu P, Lai T, Francalacci P, Curini-Galletti M. 2011b. Molecular phylogeography of the microturbellarian Monocelis lineata (Platyhelminthes: Proseriata) in the North-East Atlantic. Biological Journal of the Linnean Society 103: 117-135.

Casu M, Cossu P, Lai T, Scarpa F, Sanna D, Dedola GL, Curini-Galletti M. 2012. First evidence of self-fertilization in a marine microturbellarian (Platyhelminthes). Journal of Experimental Marine Biology and Ecology 428: 32-38.

Casu M, Scarpa F, Delogu V, Cossu P, Lai T, Sanna D, Curini-Galletti M. 2014. Biodiversity patterns in interstitial marine microturbellaria: a case study within the genus Parotoplana (Platyhelminthes: Rhabditophora) with the description of four new species. Journal of Zoological Systematics and Evolutionary Research 52: 190-202.

Cognetti G, Curini-Galletti M. 1993. Biodiversity conservation problems in the marine environment. Marine Pollution Bulletin 26 (4): 179-183.

Cohen KM, Finney SC, Gibbard PL, Fan JX. 2013 (updated). The ICS International Chronostratigraphic Chart. Episodes 36: 199-204. URL:

Cox CB, Moore PD. 2005. Biogeography: an Ecological and Evolutionary Approach, seventh edition. Oxford: Blackwell.

Curini-Galletti M, Puccinelli I. 1990. The Gyratrix hermaphro­ditus Specie Complex (Platyhelminthes: Kalyptorhynchia) in the Darwin Area (Northern Territory, Australia). Transactions of the American Microscopical Society 109: 368-379.

Curini-Galletti M, Artois T, Delogu V, De Smet WH, Fontaneto D, Jondelius U, Leasi F, Martínez A, Meyer-Wachsmuth I, Nilsson KS, Tongiorgi P, Worsaae K, Todaro MA. 2012. Patterns of Diversity in Soft-Bodied Meiofauna: Dispersal Ability and Body Size Matter. PLoS ONE 7: e33801.

Dayrat B. 2005. Toward integrative taxonomy. Biological Journal of the Linnean Society 85: 407-415.

de Queiroz K. 2005. A unified concept of species and its consequences for the future of taxonomy. Proceedings Of The California Academy Of Sciences 56, Supplement I, No. 18: 196-215.

Delogu V, Curini-Galletti M. 2009. The Parotoplana jondelii species-group (Platyhelminthes: Proseriata): a microturbellarian radiation in the Mediterranean. Contributions to Zoology 78: 99-112.

Drummond AJ, Rambaut A. 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evolutionary Biology 7: 214.

Edwards DL, Knowles LL. 2014. Species detection and individual assignment in species delimitation: can integrative data increase efficacy? Proceedings of the Royal Society B-Biological Sciences 281: 20132765.

Ezard T, Fujisawa T, Barraclough TG. 2009. SPLITS: SPecies’ LImits by Threshold Statistics. R package version 1.0-11. Available at

Fontaneto D, Flot JF, Tang CQ. 2015. Guidelines for DNA taxonomy, with a focus on the meiofauna. Marine Biodiversity 45: 433-451.

Fonseca VG, Carvalho GR, Sung W, Johnson HF, Power DM, Neill SP, Packer M, Blaxter ML, Lambshead PJD, Thomas WK, Creer S. 2010. Second-generation environmental sequencing unmasks marine metazoan biodiversity. Nature Communications 1: 98.

Fonseca VG, Carvalho GR, Nichols B, Quince C, Johnson HF, Neill SP, Lambshead JD, Thomas WK, Power DM, Creer S. 2014. Metagenetic analysis of patterns of distribution and diversity of marine meiobenthic eukaryotes. Global Ecology and Biogeography 23: 1293-1302.

Fujisawa T, Barraclough TG. 2013. Delimiting species using single-locus data and the Generalized Mixed Yule Coalescent approach: A revised method and evaluation on simulated data sets. Systematic Biology 62: 707-724.

Graff L von. 1913. Das Tierreich 35. Platyhelminthes. Turbellaria II. Rhabdocoelida. Berlin.

Jörger KM, Schrödl M. 2013. How to describe a cryptic species? Practical challenges of molecular taxonomy. Frontiers in Zoology 10: 59.

Karling TG. 1966. Marine turbellaria from the Pacific coast of North America 4. Coelogynoporidae and Monocelididae. Arkiv för Zoologi 18: 493-528.

Karling TJ. 1974. Turbellarian fauna of the Baltic proper. Identification, ecology and biogeography. Fauna Fennica 27: 1-101.

Katoh K, Toh H. 2008. Recent developments in the MAFFT multiple sequence alignment program. Briefings in Bioinformatics 9: 286-298.

Kimura M. 1980. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. Journal of Molecular Evolution 16: 111-120.

Leasi F, Norenburg JL. 2014. The Necessity of DNA taxonomy to reveal cryptic diversity and spatial distribution of meiofauna, with a focus on Nemertea. PLoS One 9: e104385.

Lemey P, Posada D. 2009. Molecular clock analysis. Pp. 362-380 in: Lemey P, Salemi M, Vandamme AM, ed., The Phylogenetic Handbook. Cambridge: Cambridge University Press.

Littlewood DTJ, Curini-Galletti M, Herniou EA. 2000. The interrelationships of Proseriata (Platyhelminthes: Seriata) tested with molecules and morphology. Molecular Phylogenetics and Evolution 16: 449-466.

Littlewood DTJ, Olson PD. 2001. Small subunit rDNA and the Platyhelminthes: Signal, noise, conflict and compromise. Pp. 262-278 in: Littlewood DTJ, Bray RA, ed., Interrelationships of the Platyhelminthes. New York: Taylor & Francis Inc.

Litvaitis MK, Curini-Galletti M, Martens PM, Kocher TD. 1996. A reappraisal of the systematics of the Monocelididae (Platyhelminthes, Proseriata) – Inferences for rDNA sequences. Molecular Phylogenetics and Evolution 6: 150-156.

Luther A. 1960. Die Turbellarien Ostfennoskandiens. I. Acoela, Catenulida, Macrostromida, Lecithoepitheliata, Prolecithophora, Proseriata. Fauna Fennica 7: 1-155.

Martens PM, Curini-Galletti M. 1994. Revision of the Archiloa genus complex with description of seven new Archilina species (Platyhelminthes, Proseriata) from the Mediterranean. Bijdragen tot de Dierkunde 64: 129-150.

Martens PM, Schockaert ER. 1986. The importance of turbellarians in the marine meiobenthos: a review. Hydrobiologia 132: 295-303.

Matthey R. 1949. Les chromosomes des Vertébrés. Lausanne: Rouge.

Meyer-Wachsmuth I, Curini-Galletti M, Jondelius U. 2014. Hyper-cryptic marine meiofauna: species complexes in nemertodermatida. PLoS One 9: e107688.

Miller MA, Pfeiffer W, Schwartz T. 2010. Creating the CIPRES Science Gateway for inference of large phylogenetic trees. Pp. 1-8 in: Proceedings of the Gateway Computing Environments Workshop (GCE). New Orleans.

Mora C, Rollo A, Tittensor DP. 2013. Comment on “Can we name Earth’s species before they go extinct?” Science 341: 237.

Norris RD, Hull PM. 2012. The temporal dimension of marine speciation. Evolutionary Ecology 26: 393-415.

Padial JM, Miralles A, De la Riva I, Vences M. 2010. The integrative future of taxonomy. Frontiers in Zoology 7: 16.

Pons J, Barraclough TG, Gomez-Zurita J, Cardoso A, Duran DP, Hazell S, Kamoun S, Sumlin WD, Vogler AP. 2006. Sequence-based species delimitation for the DNA taxonomy of undescribed insects. Systematic Biology 55: 595-609.

Posada D. 2008. jModelTest: Phylogenetic Model Averaging. Molecular Biology and Evolution 25: 1253-1256.

Puillandre N, Lambert A, Brouillet S, Achaz G. 2012. ABGD, Automatic Barcode Gap Discovery for primary species delimitation. Molecular Ecology 21: 1864-1877.

R Core Team. 2013. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0. URL:

Rambaut A, Drummond AJ. 2009. Tracer Version 1.5.0. URL:

Rosenberg NA. 2004. Distruct: a program for the graphical display of structure results. Molecular Ecology Notes 4: 137-138.

Scarpa F, Cossu P, Sanna D, Lai T, Norenburg JL, Curini-Galletti M, Casu M. 2015. An 18S and 28S-based clock calibration for marine Proseriata (Platyhelminthes). Journal of Experimental Marine Biology and Ecology 463: 22-31.

Schlick-Steiner BC, Steiner FM, Seifert B, Stauffer C, Christian E, Crozier RH. 2010. Integrative taxonomy: a multisource approach to exploring biodiversity. Annual Review of Entomology 55: 421-38.

Silvestro D, Michalak I. 2012. RAxMLGUI: a graphical front-end for RAxML. Organisms Diversity and Evolution 12: 335-337.

Sneath PHA, Sokal RR. 1973. Numerical Taxonomy: the principles and practice of Numerical Classification. San Francisco: Freeman W. H. & Son.

Sokal RR, Sneath PHA. 1963. Principles of numerical taxonomy. San Francisco: Freeman W. H. & Son

Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S. 2011. MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Molecular Biology and Evolution 28: 2731-2739.

Tang CT, Leasi F, Obertegger U, Kieneke A, Barraclough TG, Fontaneto D. 2012. The widely used small subunit 18S rDNA molecule greatly underestimates true diversity in biodiversity surveys of the meiofauna. Proceedings of the National Academy of Sciences of the U.S.A. 109: 16208-16212.

Tavaré S. 1986. Some probabilistic and statistical problems in the analysis of DNA sequences. Pp. 57-86 in: Miura RM, ed., Some mathematical questions in biology - DNA sequence analysis. Providence: American Mathematical Society.

Thomas JA, Welch JJ, Lanfear R, Bromham L. 2010. A generation time effect on the rate of molecular evolution in invertebrates. Molecular Biology and Evolution 27: 1173-1180.

Underwood AJ. 1997. Experiments in ecology. Their logical design and interpretation using analysis of variance. Cambridge: Cambridge University Press.

Vieites DR, Wollenberg KC, Andreone F, Köhler J, Glaw F, Vences M. 2009. Vast underestimation of Madagascar’s biodiversity evidenced by an integrative amphibian inventory. Proceedings of the National Academy of Sciences of the U.S.A. 106: 8267-8272.

Wagner CE, Keller I, Wittwer S, Selz OM, Mwaiko S, Greuter L, Sivasundar A, Seehausen O. 2013 Genome-wide RAD sequence data provide unprecedented resolution of species boundaries and relationships in the Lake Victoria cichlid adaptive radiation. Molecular Ecology 22: 787-798.

Will KW, Mishler BD, Wheeler QD. 2005. The perils of DNA barcoding and the need for integrative taxonomy. Systematic Biology 54: 844-851.

Zhang J, Kapli P, Pavlidis P, StamatakisA. 2013. A general species delimitation method with applications to phylogenetic placements. Bioinformatics 29: 2869-2876.

Online Supplementary Information

S1. Biometrical data (means ± SD) of populations of the Monocelis lineata complex (for details see Materials and Methods). Values are in μm. Characters: 1 = Total body length; 2 = Length of the copulatory organ; 3 = Width of the copulatory organ; 4 = Length of the penis papilla; 5 = Thickness of the proximal muscle sheath of the copulatory bulb; 6 = Thickness of the muscular components of the penis papilla; 7 = Thickness of the prostatic components of the penis papilla; 8 = Distance between vaginal pore and male pore; 9 = Distance between male pore and female pore.

S2. Results of Cochran’s C-test (C), analysis of variance (F) and Student–Newman–Keuls test (SNK) on morphological parameters of populations of the Monocelis lineata complex (for details see Materials and Methods). The significance of C-test in parameters 1 and 9, even after transformation, pointed to dishomogeneity of variances, and precluded further analysis. Characters: 1 = Total body length; 2 = Length of the copulatory organ; 3 = Width of the copulatory organ; 4 = Length of the penis papilla; 5 = Thickness of the proximal muscle sheath of the copulatory bulb; 6 = Thickness of the muscular components of the penis papilla; 7 = Thickness of the prostatic components of the penis papilla; 8 = Distance between vaginal pore and male pore; 9 = Distance between male pore and female pore.

NS = Not significant; ** = P < 0.01; *** = P < 0.001.

S3. Gaussian clustering of A) morphological, and B) karyological data, represented by barplots of individual assignment probabilities. Each individual is represented by a vertical bar, divided into coloured segments, each of which represents a cluster. Height of the segments is proportional to the probability of assignment. Barplots were constructed using Distruct 1.1 (Rosenberg, 2004).

S4. Karyometrical data (means ± SD) of populations of the Monocelis lineata complex (for details and explanation of abbreviations, see Materials and Methods).

S5. Results of Cochran’s C-test (C), analysis of variance (F) and Student–Newman–Keuls test (SNK) on karyological parameters of populations of the Monocelis lineata complex (for details and explanation of abbreviations, see Materials and Methods). a = after sqrt (x+1) transformation.

NS = Not significant; *** = P < 0.001.

S6. Reproductive outputs (offsprings × pair × week; means ± SD) of the six populations of the Monocelis lineata complex (for details and explanation of abbreviations, see Materials and Methods).

S7. Chronogram obtained by Bayesian Inference (BI) showing the interrelationships of the species based on the combined 18S + 28S D1-D6 dataset. The branch length is proportional to the number of substitutions per site. The grid indicates the mean divergence time (in myr). Node values indicate the range of divergence time (in myr).

S8. Results of species delimitation obtained by GMYC, performed on the ultrametric species tree resulting from the Bayesian analyses. The “Individuals” column shows the number of specimens for each found entity. The “Taxon ID” column shows how many specimens have been assigned to the corresponding entity.

S9. Results of species delimitation obtained by PTP, and its Bayesian implementation bPTP, performed on the Bayesian phylogenetic species tree. The “Individuals” column shows the number of specimens for each found entity. The “Taxon ID” column shows how many specimens have been assigned to the corresponding entity.

S10. Results of species delimitation obtained by ABGD on the combined dataset (18S + 28S D1-D6) using K2P genetic distances (Kimura, 1980). The “Individuals” column shows the number of specimens for each found entity. The “Taxon ID” column shows how many specimens have been assigned to the corresponding entity.

S11. Results of K/θ applied to nine pairs of sister clades (see Fig. 1( for topology details), representative of the taxa here discussed (i.e., Monocelis lineata complex sensu stricto, M. exquisita, M. algicola, M. fusca). The “Taxon ID” column shows clades which were compared. In order to be as conservative as possible, for each ratio the upper values of θ have been used. Therefore, asterisks in the “K/θ” column indicate which values have been considered for the choice of the results. Clades that showed values of K/θ ratio has ≥ 4 are considered as entities with a 95% probability of having an independent evolutionary history.

S12. Map of the Porto Puddu area (north Sardinia) (pointed by the arrow in the inset), showing sampling sites and position of the pier (arrow) at the entrance of the lagoon. : present known distribution of Monocelis exquisita nov. sp.; *: sites where the species has not been found since 2004; O: species absent in 2004 and 2015.



Appendix 1

List of populations sampled. Accession numbers of sequences refer to GenBank codes. G = number of specimens used for genetics analyses. a The suffixes “o” and “x” indicate the presence and the absence of the eye pigment, respectively. b The codes KFo and TJx here substitute ICo and SVx in Casu and Curini-Galletti (2004), respectively. c Sibling codes used in Casu and Curini-Galletti (2004). d Abbreviation “Br” and “Ma” indicate brackish-water and marine populations, respectively. e Superscript codes refer to the specific analysis applied to individuals: M = morphology; K = karyology; R = reproductive biology; G = genetics. Accession numbers of newly sequenced taxa are in bold.


Appendix 2

In the following descriptions, Pore Index (a/c: b/c: c) is based on the following distances: a = mouth - vagina; b = vagina - male pore; c = male pore - female pore, with c = 1. Karyotype formula is represented as follows: fundamental number (NF) (according to Matthey, 1949); relative length and centromeric index of each chromosome; chromosome nomenclature within parentheses (m = metacentric; sm = submetacentric; st = subtelocentric; t = acrocentric).

Molecular diagnostic characters were checked on a subset composed by specimens belonging to the two new species, M. lineata complex sensu stricto, and M. fusca, in which molecular pure characters have been underlined (see Tables 1, 2). The state of the pure characters is present only across all members of a single clade (see Jörger and Schrödl, 2013).

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 The LSID for this publication is: The electronic edition of this work was published in a journal with an ISSN, and has been archived.

Types have been deposited at the Swedish Museum of Natural History (SMNH) (Stockholm, Sweden) and in the collections of the Zoological Museum (CZM), University of Sassary (Italy).

Phylum Platyhelminthes Minot, 1876

Order Proseriata Meixner, 1938

Family Monocelididae Hofsten, 1907

Subfamily Monocelidinae Midelburg, 1908

Genus Monocelis Ehrenberg, 1831

Monocelis algicola Curini-Galletti and Casu nov. sp. EC8B9EF7C0AF

(Fig. 5A-D()


Fig. 5A-D: Monocelis algicola nov. sp.: sagittal reconstruction of the genital area (A); general organisation of a live specimens (B); sagittal section of copulatory organ (C); spermatogonial metaphase plate; arrows point to heterobrachial Chromosome II (D). E-H: Monocelis exquisita nov. sp.: sagittal reconstruction of the genital area (E); general organisation of a live specimens (F); sagittal section of copulatory organ (G); spermatogonial metaphase plate; arrows point to Chromosome III (H). I-M: Monocelis lineata s.s.: sagittal reconstruction of the genital area (I); sagittal section of copulatory organ (H); spermatogonial metaphase plate (M), from specimens from Helsingor, Denmark (HEo). Scale bar : C, G, L = 10 μm; D, H, M = 5 μm. Abbreviations used in figures: b: bursa, co: copulatory organ, e: eye, fd: female duct, fg: female glands, fp: female pore, gl: gut lumen, ml: muscular lining, mp: male pore, ov: ovary, pg: prostatic glands, ph: pharynx, pp: penis papilla, pt: prostatic tissue, rh: rhabdoid gland, sp: sphincter, st: statocist, t: testis, v: vagina, vi: vitellaria, vs: seminal vesicle.

Holotype. Italy: Sardinia, Cala Rossa (Lat. 41.027595N; Long. 8.892886E), 2-3 m deep on algae (05.15.2002): one specimen sagittally sectioned (SMNH-Type-8770).

Paratypes: same data as holotype, 32 specimens sagittally sectioned (CZM 614-645).

Other material: 15 specimens studied karyologically (in non-permanent mounts); five specimens sequenced.

Etymology: the specific epithet refers to the habitat of the species (from latin alga: seaweed; colere: to dwell).

Description. A very slender Monocelis, slightly tapering anteriorly. Living specimens very active; extremely adhesive due to the presence of numerous adhesive glands in the tail. The fixed holotype is 2.8 mm long: living, extended specimens up to 3 mm long. Many specimens show faint, irregular brownish lines in the cephalic area. With a large, pigmented eye-shield, just in front of the statocist. Epithelium with insunk nuclei, ciliated (cilia about 4 μm long) apart from the dorsal side of the caudal tip. With two types of rhabdoid glands: a) large, up to 30 μm long, numerous dorsally and caudally, with fine, filamentous content, appearing dot-like in transverse sections; b) smaller (up to 15 μm), more globular glands, with much thicker content (about 0.5 μm in diameter), particularly abundant caudally. With well developed ventral longitudinal musculature; rest of subepidermal musculature poorly developed. With a short, tubular pharynx, about 1/9th total body length, almost at mid-body. The pharynx is ciliated (cilia 3-4 μm long) a part from distal tip, where pharyngeal glands discharge. Musculature layers inverted with respects to body musculature, poorly developed a part from the layers surrounding the inner lumen. Oesophagous 1/4th to 1/5th the total length of the pharynx.

Male genital system. With about 30 testes, irregularly arranged in front of the pharynx. The copulatory organ is spherical (48.8 ± 1.03 μm wide, 47.18 ± 0.93 μm high; 30 measurements) (Fig. 5C(). It consists of a seminal vesicle, surrounded by a thin muscular layer, and provided with a short penis papilla. Thin (up to 2 μm), proximal musculature mostly longitudinal. Penis papilla 7.31 ± 2.01 μm long, surrounded by a thin muscular coating. The proximal basis of the papilla is lined by stronger, circular muscles. The penis papilla is lined internally by a thick layer of prostatic tissue, whose cell nuclei are mostly located outside the copulatory organ. Prostatic secretion dot-like. Male atrium small; opening to the outside though the male pore.

Female genital system. Ovaria just in front of pharynx. Vitellaria run from the level of the first testes to in front of bursa. Oviducts join behind pharynx to form an unciliated female duct. Just in front of the copulatory organ, the female duct is connected ventrally to a short, unciliated vaginal duct (10-15 μm long), surrounded by a strong sphincter, and opening to the outside through a vagina. Immediately above this duct, a very short bursal duct, also surrounded by a strong muscular sphincter, leads to a large bursa, up to 100 μm across, consisting mostly of resorbing vacuoles, many of which abutting the gut, with sperm in different stages of degeneration. Distally to this vaginal-bursa complex, the female duct is lined with an irregular, high, at least partly resorbiens epithelium, runs over the copulatory bulb and opens to the outside just posterior to the male pore through a female pore, surrounded by female glands (Fig. 5A().

Karyotype. Chromosome number n=3. With a distinctly heterobrachial chromosome pair. Chromosome I and II nearly even in size; Chromosome III about 83% the size of Chromosome I. Karyotype formula: FN=5; Chromosome I: 35.69 ± 0.56; 42.19 ± 0.8 (m); Chromosome II: 34.88 ± 0.62; 9.26 ± 0.78 (t); Chromosome III: 29.46 ± 0.43; 44.15 ± 0.61 (m) (Fig. 5D().

Diagnosis. slender Monocelis species with a large, pigmented eye-shield. With large rhabdoids with thin fibrillar content, and smaller rhabdoids with coarser content. With about 30 testes. Spherical copulatory organ (about 48 μm across), with a thin muscular coating, and a poorly developed penis papilla, with a thick layer of prostatic tissue. Prostatic secretion finely granular. With a large, mostly resorbiens bursa, and a vagina close to male pore; female duct at least partly resorbiens distally. Karyotype with n=3, and Chromosome II acrocentric. Pore Index: 6.7: 0.9: 1.

Molecular diagnosis: see Table 1.

Monocelis exquisita Curini-Galletti and Casu nov. sp. 895285647043

(Fig. 5E-H()

Holotype. Italy: Sardinia, Porto Puddu (Lat. 41.192306N; Long. 9.340800E), lagoon, about 30 cm deep in coarse silty sand (07.01.2002): one specimen sagittally sectioned (SMNH-Type-8769).

Paratypes: same data as holotype, 29 specimens sagittally sectioned (CZM 646-674).

Other material: 15 specimens from type locality studied karyologically (in non-permanent mounts); five specimens sequenced from type locality. Porto Puddu, mouth of creek close to the lagoon (see map in Fig. S3) (Lat. 41,983; Long. 9,317), about 20 cm coarse silty sand (07.01.2002; 05.28.2014), numerous specimens studied alive, three specimens studied karyologically.

Etymology: the specific epithet reflects the extensive search for populations of the species outside the type locality (lat. exquisitus: searched for, sought out).

Description. A slender, active Monocelis. The fixed holotype is 2.5 mm long: living specimens up to 3 mm long. Some specimens show faint, irregular brown lines in the cephalic area. With a pigmented eye-shield in front of statocist, irregular in shape; in a few specimens split into two closely lying spots. Epithelium with insunk nuclei, ciliated (cilia about 3-3.5 μm long) apart from the dorsal side of the caudal tip. With large rhabdoid glands, particularly numerous dorsally and caudally, in some areas appearing as a continuous layer parallel to the surface, with coarse-grained, thick (> 1 μm in diameter) fibrillar content, intensely stained with hematoxylin (Fig. 5C(). A few, scattered, subepithelial glands, with fine-grained, dot-like eosinophilous content are present dorsally. Posterior end with adhesive glands. Pharynx at the beginning of the posterior half of body, comparatively long (between 1/7th and 1/8th of total body length). It is ciliated externally (cilia 2-2.5 μm long), and internally (cilia 3-4 μm long) in its distal half, a part from its tip, where pharyngeal glands discharge. With very strong inner circular musculature. Oesophagous about 1/4th the total length of pharynx.

Male genital system. With numerous (40-60) testes, irregularly arranged in front of the pharynx. The copulatory organ is nearly spherical (44.12 ± 0.96 μm wide, 39.02 ± 0.92 μm high; 30 measurements), backward oriented, with an angle of about 40° to the vertical axis (Fig. 5G(). It consists of a seminal vesicle, surrounded by a thick muscular layer, with a thin, pointed penis papilla, ranging 12.86 ± 0.43 μm in length (30 measurements). Proximally, the muscular coating, up to 3 μm thick, consists mostly of longitudinal muscles. Circular musculature more developed distally; penis papilla lined by a thick layer of outer circular muscles, with a thin layer of longitudinal, inner musculature. The distal part of the copulatory bulb is lined internally by a thin layer of glandular epithelium, whose cell bodies are placed mostly externally to the copulatory bulb. These glandular cells have a fibrillar content. Male atrium small.

Female genital system. Similar to the previous species in general arrangement. Ovaria just in front of pharynx. Vitellaria run from the level of the first testes to in front of bursa. The vaginal duct is short (about 15 μm long), unciliated and surrounded by a strong sphincter. Above it, a large vacuolar bursa (up to 150 μm in diameter), abutting the gut, is present. The female duct, lined by an irregular, high, at least partly resorbiens epithelium, opens to the outside just posterior to the male pore through a female pore, surrounded by female glands (Fig. 5E().

Karyotype. Chromosome number n=3. Karyotype with Chromosome I and II almost of the same length; Chromosome III markedly smaller, about 60 % the size of Chromosome I. Karyotype formula: FN=6; Chromosome I: 39.26 ± 0.51; 45.07 ± 0.7 (m); Chromosome II: 37.38 ± 0.46; 31.72 ± 1.6 (sm); Chromosome III: 23.35 ± 0.42; 29.43 ± 1.09 (sm) (Fig. 5H().

Diagnosis. Slender Monocelis species with pigmented eye-shield. With large rhabdoid glands with coarse fibrillar content, and smaller rhabdoids with dot-like, eosinophylous content. With 40-60 testes. With a nearly spherical, backward oriented copulatory organ (about 40 μm across), with a slender penis papilla about 13 μm long, internally coated by a thin layer of prostatic tissue. Prostatic secretion fibrillar. With a large, mostly resorbiens bursa and external vagina close to male pore. Karyotype with n=3, and Chromosome III markedly smaller than the other pairs. Pore index: 6.3: 1.1: 1.

Molecular diagnosis: see Table 2.

Species justification

In addition to diagnostic characters of sequences (see above), the two new species are uniquely characterised by Chromosome III with low value of centromeric index and markedly smaller than the other pairs, and the penis papilla with a very thin prostatic component (Monocelis exquisita nov. sp.); Chromosome II acrocentric, extremely thin layer of proximal musculatura of the copulatory bulb, and the penis papilla with a thick prostatic component (M. algicola nov. sp.) (Tables S1-S4). The two new species appear externally quite similar, as both are more slender than the rest of the species-group, and may present similar cephalic pigmentation. In addition to the remarkably different morphology of their copulatory bulbs, the two species differ for morphology of prostatic secretion (dot-like in M. algicola nov. sp.; fibrillar in M. exquisita nov. sp.), and content of the larger, dorsal, rhabdoid glands: thin, fibrillar in M. algicola nov. sp., thick, coarse-grained in M. exquisita nov. sp. It is worth mentioning that the morphology of both rhabdoid and prostate gland secretion found in M. exquisita is also found in HEo (close to type locality of Monocelis lineata, see above) (Fig. 5I() and is then possibly plesiomorphic.