Morphological datanext section
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.
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).
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.