We documented the main taxonomic characters of the tubular shells (Fig. 1), such as the morphology of aperture, septum, and mucro, and measured size and diameter of the shells. Initial species identification in the field was carefully revised in the laboratory, and specimens were assigned to species according to these shell characteristics. The microsculpture of the shell of one representative of each putative morphospecies was investigated via scanning electron microscopy (SEM), whenever a voucher was available (Table 1).

Microscopic debris on the shell was manually removed using an eyelash, and the shell rinsed in 96% ethanol. Specimens were dried by evaporation of the ethanol and transferred onto SEM stubs covered with self-adhesive carbon stickers. We used a sputter coater Polaron SC510 to coat the samples with gold in argon atmosphere. The shells were analyzed with a LEO 1430 VP SEM at a voltage of 15 kV.

All light microscopic images and SEM-micrographs are available through FigShare (https://figshare.com/projects/Central_American_Caecidae/84929).

DNA was extracted from 121 of the 132 investigated specimens. The sputter-coated individuals previously investigated by SEM were crushed mechanically using pestles (Bergmeier et al. 2016); specimens investigated only by LM were also crushed if tissue was not already separated. Subsequently, DNA was extracted by the procedure of Knebelsberger and Stöger (2012) combining lysis with 2-mercaptoethanol in CTAB buffer, chloroform-isoamyl precipitation, and recovery using columns with silica-membrane (Nucleo Spin, Macherey-Nagel GmbH & Co. KG, Düren, Germany). The DNA was eluted twice with 25 µl aliquots of pre-heated elution buffer to gain high yield. DNA of specimens deposited to the USMN-Smithsonian Institution were extracted in the Laboratories for Analytical Biology, SI using the standard protocols of the Autogen Prep 956 Extractor (eluting with 100 µl Autogen R9 buffer). Three different markers were partially amplified by PCR: mitochondrial cytochrome oxidase subunit I (COI) and 16S rRNA gene, and nuclear 28S rRNA gene, using the standard PCR primers for gastropods (see Klussmann-Kolb et al. 2008). We either used the Phire polymerase (Thermo Fisher Scientific Inc., Waltham, USA) with the following protocol for PCRs on 16S/ COI resp. 28S at the LMU: 98 °C 1 min (98 °C 30 sec; 46–50 °C 20 sec; 72 °C 20 sec) × 36–38 cycles 72 °C 1 min resp.98 °C 30 sec (98 °C 15 sec; 55–60 °C 5 sec; 72 °C 20 sec) × 35 cycles 72 °C 1 min or the KlenTaq polymerase (AB Peptides, Inc.) with the following program for sequences generated at the SI: 95 °C 3 min, (95 °C 30–45 sec; 48–52 °C 30–45 sec; 72 °C 45–90 sec) × 35–40, 72 °C 7 min. PCR products were either cleaned using a spin column purification kit (Zymo Research, Irvine, California, USA) or were purified with QIAquick (Qiagen Inc.). Samples at the LMU were cycle sequenced on an ABI 3730 48 capillary sequencer (Applied Biosystems, Foster City, USA) using Big Dye 3.1 (Thermo Fisher Scientific Inc.) at the sequencing service of the Ludwig-Maximilians-Universität (LMU) Biocenter, Munich, Germany. At the SI, cycle sequencing was also conducted with BigDye chemistry (PerkinElmer) and standard cycles (4 min denaturation at 96 °C, followed by 25 cycles of 10 sec at 96 °C, 5 sec at 50 °C and 4 min at 60 °C), and sequenced on an ABI 3730xl 96-well capillary sequencer. In total, 34%, 43% and 50% of the partial COI, 16S rRNA, and 28S rRNA gene sequences, respectively, were successfully amplified and sequenced. All sequences were edited in Geneious Prime (vers. 11.02011, Biomatters, Ltd., Auckland, New Zealand). Primer sequences were removed and base calls checked for misreads against their chromatogram. The sequences were then compared to sequences in the public database NCBI GenBank (http://ncbi.nlm.nih.gov/genbank) by using the BLAST online web service to check for putative contamination. In total 29 COI, 40 16S rRNA and 43 28S rRNA gene sequences were deposited in NCBI GenBank (see Table 1 for accession and voucher numbers).

Multiple sequence alignments of the 28S rRNA and COI genes were constructed using Mafft (vers. 7.419; Katoh et al. 2002; Nakamura et al. 2018) with default parameter settings. The mitochondrial 16S rRNA sequences were aligned using the program Muscle (vers. 3.8.31; Edgar 2004) with default parameter settings. Alignments were visualized using Seaview (vers. 3.2; Gouy et al. 2009). COI sequences were translated into amino acids. The program Gblocks (vers. 0.91b; Castresana 2000; Talavera and Castresana 2007) was applied to the 16S and 28S rRNA gene alignments to check for unambiguously aligned sites. Proposed exclusion sites were reviewed, adjusted, and subsequently removed (alignments before and after editing are deposited at https://doi.org/10.5281/zenodo.3613958). Sequences available from NCBI GenBank for in-group taxa (C. glabrum (Montagu, 1803) and C. glabellum (A. Adams, 1868)), as well as for out-group taxa were added (Table 3). Outgroups were assigned based on the recent molecular phylogenetic analyses by Golding (2014a, b) and Criscione and Ponder (2013).

Two combined data sets were generated: (1) a concatenated alignment of all three marker genes and (2) a concatenated alignment comprising only the mitochondrial 16S rRNA gene, and COI. The data were combined into single matrices using P4 (Foster 2004). The combined data sets were then partitioned by gene and COI codon position.

Maximum likelihood (ML) and Bayesian inference (BI) were used to construct the phylogenetic tree from single genes and from combined and partitioned alignments. For each alignment jModelTest2 (vers. 2.1.10; Darriba et al. 2012) was run and the calculated likelihood scores weighted under the Akaike Information criterion (AICc) (Hurvich and Tsai 1989) which suggested GTR+I+G as the best fitting model. ML was performed using IQ-TREE (multicore vers. 1.6.7.1 for Linux 64-bit; Nguyen et al. 2014) with the GTR+G4+FO model (equivalent to GTR+G in RAxML vers. 8.2; Stamatakis 2014) with 300 bootstrap replicates. Bayesian MCMC analyses were performed using the program MrBayes (vers. v.3.2.6; Huelsenbeck and Ronquist 2001) with the same model. The Bayesian analyses were run in duplicates by default, with each run having four parallel Markov chains (MCMC) to estimate posterior probability support. Each chain was run for 5 million generations, sampling trees every 1000^th generation. Sampled trees were combined into a consensus tree after the first 1000 sampled trees (1000000 generations), considered as ‘burn-in’, were discarded. A general time-reversible model of nucleotide substitutions with a gamma-distribution of among-site rates (GTR+G) was used for the ML analyses. All trees were visualized and annotated using Figtree (vers. 1.4.4; Rambaut 2007). Boostrap support values (BS) > 85% and posterior probabilities > 0.95 were considered statistically significant.

Four different methods of species delineation were used with both the COI and 16S rRNA gene mitochondrial data sets. The Automatic Barcode Gap Discovery (ABGD) webserver was used to partition the data set into putative species based on the calculated gap between intra- and interspecific genetic differences (https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html; Puillandre et al. 2012). J ModelTest2 (vers. 2.1.10; Darriba et al. 2012) was applied to the uncorrected COI and 16S rRNA gene alignments and the parameters were weighted under the corrected Akaike Information criterion (AICc) (Hurvich and Tsai 1989). For both alignments the Jukes-Cantor (JC69) as well as Kimura (K80) model showed to be within the 100% confidence interval however, K80 had slightly higher likelihood scores. Both models were applied with the default settings (TS/TV = 2.0, relative gap width = 1.5, Pmin = 0.001, and Pmax = 0.10) on the uncorrected COI and 16S rRNA gene alignments.

To evaluate haplotype connectivity, we generated haplotype networks based on the COI as well as the 16S rRNA gene sequence alignment using the software TCS (vers. 1.21; Clement et al. 2000) using the standard 95% parsimony setting. Ambiguous sites in both alignments were removed to prevent the creation of artificial haplotypes.

The bPTP web server (https://species.h-its.org/) was used to conduct the Bayesian implementation of the PTP model for species delimitation (Zhang et al. 2013) on the optimal ML trees of the individual and the combined COI and 16S rDNA datasets. We applied the default settings with 100000 generations, thinning for each 100^th sample with a burn-in of 10% and checked for convergence of the MCMC chains of each run. Posterior probability (PP) support values above 0.95 were considered as strong support.

For the General Mixed Yule-Coalescent model (GMYC) (Pons et al. 2006), ultrametric trees from the COI, 16S rRNA gene, and combined COI and 16S rRNA gene data were obtained using a time calibrated Bayesian evolutionary analysis in Beast (vers. 1.7.4; Drummond and Rambaut 2007). For the tree prior, we used a Yule process and two fossil records, Caecum cooperi and Caecum imbricatum [2.58–1.80 myr] (Mansfield 1930; Cooke 1936; Ward and Blackwelder 1987) and the in-group Caecum [50–55 myr] (Goedert and Raines 2016) with a lognormal distribution (logL). The analysis was run with the GTR substitution model and under a strict clock assumption. The analysis was started from a random tree and two Markov chains run for 10 000 000 generations with a sampling frequency of 1000. Convergence of the chains was checked in Tracer (vers. 1.7.4.; Rambaut et al. 2018) and effective sampling sizes (ESS) were confirmed as > 200 for all values (Rambaut et al. 2018). The first 10% of sampled trees were removed as burn-in and the trees were combined in TreeAnnotator (vers. 1.7.4.; Drummond et al. 2012) using the maximum clade credibility option and mean node height. The ultrametric trees were uploaded to the web server (https://species.h-its.org/gmyc) for single, as well as, multiple threshold GMYC analyses.

The software QUIDDICH (vers. 1.0.0; Kühn and Haase 2019) was used to identify the diagnostic molecular characters of morphologically cryptic species. We extracted diagnostic characters of type 1 (i.e., characters, which distinguish each individual of the investigated species from other caecids with a fixed character state in the investigated species) and of type 2 (i.e., characters, which distinguish each individual of the investigated species from all other caecids, but vary also within the investigated species) from the COI, 16S rRNA gene and 28S rRNA gene alignments of the same dataset also used for the species delineation and phylogenetic analyses.

In our phylogenetic analyses Caecidae form a well-supported clade (1.0 PP, 99% BS; Fig. 3). The two established genera Caecum and Meioceras, however, are not recovered as reciprocally monophyletic but instead species of Meioceras group among Caecum species in different parts of the tree (Fig. 3, taxa highlighted in yellow): M. nitidum sister to C. heptagonum (0.96 PP), and M. cubitatum sister to C. cf. semilaeve (no statistical support). The phylogeny groups the Caecidae into 21 clades which show moderate to high support values ranging from 0.95 PP/85% BS to full support (Fig. 3). Other clades are only statistically supported by one analysis (C. regulare, 92% BS) or do not have statistical support (C. donmoorei). The sister group relationships of C. pulchellum and C. regulare (1.0 PP, 95% BS), and C. cooperi and C. imbricatum (1.0 PP, 100% BS) are well supported; otherwise, deeper nodes and higher-level relationships among clades are not supported. In agreement with the molecular data, C. pulchellum and C. regulare as well as C. cooperi and C. imbricatum show morphological similarities in shell ornamentation and microsculpture. Inconspicuous specimens with smooth shells and few characters that were morphologically ascribable to C. glabrum or the American Pacific look-alikes like C. glabriforme are polyphyletic and form four lineages separated by branches of comparable length to morphologically distinct species (Fig. 3, highlighted in blue). These lineages are distinct from C. glabrum from the North Atlantic (Table 3) included in the analyses, indicating the presence of morphologically cryptic species in this ‘C. glabrum species complex’.

Figure 3. 

Optimal ML tree of the concatenated 28S rRNA, 16S rRNA and COI genes partitioned by genes and COI codon positions. Bootstrap values (below nodes) of the ML analysis are shown for values > 80% and posterior probability support (above nodes) of the BI analysis are shown for values > 0.95. Specimens previously classified as Meioceras are indicated in yellow color. Smooth, translucent specimens, lacking diagnostic features and summarized in the ´Caecum glabrum-like complex` are indicated in blue. C. = Caecum, M. = Meioceras, MOTU I/ MOTU II = molecular operational taxonomic unit within the ´Caecum glabrum-like complex`. Figured specimens are all to the same scale. Scale bar: 1 mm.

The methods that were used for species delineation are largely congruent with regard to the assignment of taxa to molecular operational taxonomic units (MOTUs), however individual analyses deviate and evidently differences occur due to incomplete sampling of one of the markers (Fig. 4, Table 1). Both PTP/ bPTP and GMYC (single threshold) delimit 21 MOTUs for the concatenated dataset of COI and 16S rRNA genes (excluding the species whose sequences were retrieved from NCBI GenBank, i.e. North Atlantic C. glabrum and Japanese C. glabellum). These results are in concordance with the preliminary species hypotheses based on morphological investigation and the molecular phylogenetic tree (Fig. 3) with the exception of additional splits of C. debile and C. regulare into two distinct MOTUs each. Caecum cf. teres resulted in a single species for 16S rDNA alone, and the Bayesian implementation of bPTP split M. nitidum into two separate species based on the 16S rRNA genes (however, support value for the split is 0.501%). The multiple threshold analyses in GMYC additionally splits C. invisibile sp. nov. of the ‘C. glabrum complex’ into two MOTUs, as does TCS but into differing entities. In analyses of individual datasets (numbers not directly comparable due to missing data) ABGD identified 15 MOTUs in our 16S rRNA gene dataset (Fig. 4), while the COI dataset resulted in a hypothesis of 10 MOTUs independent of the application of the JC69 or the K80 model. In comparison to the other methods, TCS appears to oversplit MOTUs (see e.g., TCS analyses of the COI of C. donmoorei in Fig. 4). The algorithm of this haplotype-network software splits the 16S rDNA dataset into 19 independent haplotype networks, while it recovered 13 networks for the COI dataset (Fig. 4). Additionally, TCS also splits MOTU II of the ‘C. glabrum-like complex’ into two networks and C. donmoorei into three independent networks based on 16S rRNA sequence data. Haplotype networks divided C. cf teres and C. cf. strangulatum into two unconnected networks. However, the split is not congruent with the two monophyletic sister populations of the species tree (Fig. 4). In summary, we consider only splits relevant, which are supported by at least two different analyses or markers, singular deviating signal might either resemble errors in analyses or might be informative in population analyses (for more details see remarks in Systematics section).

Figure 4. 

Molecular based species delimitation of Central American Caecidae. Guide tree used for PTP and bPTP based on the optimal likelihood tree of the concatenated three-marker dataset. Color codes indicate our preliminary species hypothesis derived from the phylogenetic tree. Color bars reflect the species delimitation suggested by the four consulted species delimitation programs (including ML and Bayesian implementation for PTP and single and multiple threshold for GMYC). Bars are missing where no sequence data obtained.

Twenty-four specimens from Central American waters are smooth and glossy without ornamentation except for occasional growth lines (i.e., possess few shell characteristics), but vary in adult shell length between 0.7 and 2.5 mm (Figs 3, 5F–L, O–R). Morphologically, these specimens all closely resemble Caecum glabrum (Montagu, 1803) which is one of the best-known species of caecids, and abundant in the northern Atlantic (Montagu 1803; Wood and Harmer 1848; Götze 1938; Chambers 2009). Caecum glabrum was originally described from Biddlesford Bay and Barnstable, Devon, England, the included sequences from GenBank (see Table 3) originates from specimens collected in Norway, but own unpublished data from Roscoff, northern France, supports the wide distribution range of C. glabrum along European coastlines based on molecular data. We refer to cryptic species with simple shells lacking characteristic features as ‘Caecum glabrum-like’ species complex. In previous works, specimens similar to C. glabrum have also been described from the Pacific (C. glabellum as Brochina glabella A. Adams, 1868 from Akashi, Japan, C. glabriforme Carpenter, 1857 and C. corrugulatum Carpenter, 1857, both from Mazatlán, Mexico and C. parvulum de Folin, 1867 from Panama Bay, Brazil), have been reported and described from Japan, Hawaii, and central America (Carpenter 1855–1857; Adams 1868; Lightfoot 1993b; Pizzini et al. 2007, Takano and Kano 2014). Our study clearly shows an independent evolutionary origin of C. glabrum from the northeast Atlantic and C. glabellum from Japan, and the cryptic C. glabrum-like MOTUs from central America (see Fig. 3, ‘C. glabrum-like-complex’ highlighted in blue color). Our molecular species delineation revealed a minimum of four cryptic MOTUs (see above, Fig. 4).

The Caecidae are currently classified in ten genera (MolluscaBase 2019) of which two, Caecum and Meioceras, can be found in the Central American region. We investigated two species classified as Meioceras and 15 Caecum species, including one species new to science (C. invisibile sp. nov.), and two candidate species (MOTU I and II). Three individuals that were originally assigned to Meioceras are resolved among species of the genus Caecum in our molecular phylogenetic analyses and, moreover, have independent evolutionary origins within Caecum (Fig. 3). Our data confirms the existence of at least two valid Meioceras species (i.e., M. nitidum and M. cubitatum), which however, based on our data should be reassigned genus Caecum. For taxonomic stability, we refrain to reallocate these species at present, until the molecular sampling can be expanded and further data is available supporting our initial results. Unfortunately, we lack material of M. cornucopiae (the type species by subsequent designation) and of a putative fourth species, M. tumidissimum, both described from Brazil (de Folin and Périer 1869). These taxa are needed to settle the debate on the number of species and to clarify the validity of the genus Meioceras. Our findings indicate, however, that the more bulbous shell in Meioceras when compared to a more tube-like shell in Caecum might not justify generic subdivision. Moreover, the diagnostic spiral growth pattern in the larval and juvenile shell of Meioceras (Bandel 1996) might not be used for unambiguous discrimination as our phylogeny indicates that such patterns have evolved independently at least two times within the genus Caecum. The bulbous adult shell with an oblique constricted aperture was thought to develop from the preceding ontogenesis of a helicoidal shell section, observable at the beginning of the second growth stages (De Folin 1880; Absalão and Pizzini 2002). This diagnosis has been problematic because in the past other species (Absalão and Pizzini 2002) that also express a similar shell-shape had been classified as Caecum due to a lack of observations of the juvenile stadia or the lack of the aforementioned growth pattern (e.g., Caecum ryssotitum, Bandel 1996). By contrast, species with typical tube-shaped Caecum-shells are known to also have curved growth axes (e.g., C. antillarum and C. japonicum referred to as C. glabellum in Bandel 1996: 65, figs 8, 9). Differences between these growth patterns and those of Meioceras seem to be negligible. Thus, it remains to be tested whether the remaining two Meioceras species form a monophylum separate form Caecum and whether alternative diagnostic morphological or molecular characters can be found to justify the generic subdivision of Western Atlantic Caecidae. Alternatively, such studies may confirm that Meioceras is a junior synonym of Caecum. The present study might also have consequences for the recently established genus Mauroceras, which unites Indo-Pacific Caecidae formerly classified as Meioceras (Vannozzi 2019). But, in contrast to Caecum and Western Atlantic Meioceras, which cannot be clearly separated based on variable growth patterns, Mauroceras is diagnosed by a planorbid protoconch with a clear sinusigera, which at present justifies its generic status.

The taxonomy of Central American Caecidae has been based on macroscopic shell characters and, consequently, type-species are often poorly defined, and has made the established taxonomy prone to multiple descriptions of synonyms and the establishment of ambiguous species-complexes that are typical for many clades of micromolluscs (Golding 2014b). Modern microsculptural analyses have greatly increased the reliability of shell-based taxonomy and the availability of diagnostic characters in the otherwise largely featureless caecid shells (Pizzini et al. 2013; Vannozzi 2017). However, distinct shells based on coarser diagnostic features can have a similar microsculpture (Vannozzi 2017), suggesting that shell microsculpture should be co-evaluated with traditional diagnostic features and, indeed, that it might be especially valuable to discriminate closely related species. In Central American Caecidae the presence of a series of morphologically highly-similar ribbed taxa (i.e., C. compactum Carpenter, 1857, C. quadratum, C. clathratum, C. gurgulio Carpenter, 1858, C. pulchellum, C. regulare, Caecum planum de Folin, 1874) with controversial species status and inconsistent synonymization (Moore 1972) are especially problematic. Here we report SEM-based shell microsculpture that can distinguish taxa and justify the independent species status of C. regulare and C. pulchellum (compare fusiform lobes (Fig. 7C) with fine longitudinal lamellae (Fig. 7A). However, C. clathratum does not possess a unique microsculpture and we are lacking SEM data for our investigated specimen of C. donmoorei. Nevertheless, molecular species delineation analyses confirm the existence of four genetically distinct species among those Central American ribbed caecids (i.e., C. pulchellum, C. regulare, C. donmoorei and C. cf. clathratum), highlighting the value of complementary molecular analyses to detect possible synonyms or confirm the validity of existing species in taxonomically problematic species complexes.

The different growth stages of caecid development present an additional problem for taxonomic circumscription, which cannot be overcome easily by microsculptural analyses because, the shape and some patterns of ornamentation appear late in development. This often results in the incorrect assignment of different growth stages even at the generic level (Absalão and Pizzini 2002). In consequence, it requires time consuming comparisons of hundreds of shells for reliable species description and identification (Lightfoot 1992a, 1993a), unfeasible in modern times of taxonomic impediment. The molecular analyses presented in this study show that barcoding markers (i.e., partial mitochondrial COI and 16S rRNA genes) are a valuable tool to address the challenges of caecid taxonomy and that molecular species delineation analyses can reliably identify groups of closely related specimens, therewith providing objective data on intraspecific variability of shell characters. Above all, they enable an unambiguous assignment of juvenile forms in different growth stages to their fully developed adult morphologies (see e.g., C. heptagonum in Fig. 5G, C. pulchellum and C. cf. teres in Fig. 7A, D). Based on a purely morphological approach, these juveniles would have remained unidentified and unaccounted for in biodiversity data, and their contribution to caecid diversity would have been lost. However, in some cases, juveniles could not be matched to their adult counterparts using molecular data since we had no adult animals in our sample. These taxa identified by the molecular data could not be named (e.g., Caecum sp. MOTU II). These examples highlight, how the successful identification of juveniles lacking morphological diagnostic features by means of their genetic fingerprints requires an extensive barcode library of Central American Caecidae as a taxonomic reference. The barcodes of the morphospecies investigated here are the first contribution to such a reference library that can help to provide a baseline and enhance future identification. In general, the poor taxonomic coverage of gastropods and marine invertebrates in public molecular databases such as NCBI GenBank has been identified as a major obstacle to making effective use of molecular barcoding approaches (e.g., to assign spawn to adult specimens; Puillandre et al. 2009). Thus, it is hoped that in the future the scientific community will be able to invest more of its financial and personnel capacities in integrative faunistic approaches that strengthen fundamental biodiversity research.

In biodiversity assessment and conservation biology, molecular species delineation has also demonstrated its potential for identifying cryptic species (Bickford et al. 2007; Jörger et al. 2012; Lemer et al. 2014; Leasi et al. 2016). In revealing cryptic taxa, our study indicates that the species diversity of caecids may have been underestimated until now. Unsurprisingly, the cryptic species, which we identified, are those of particularly small, feature-poor, caecids with few diagnostic characters (see Fig. 5C–F). Indeed, our analyses suggest that meiofaunal character-poor caecids (assigned to the ‘Caecum glabrum-like’ species complex) have evolved several times independently from the larger ornamented caecids in the Central American region. The same may have happened in the northern Atlantic C. glabrum and Northwest Pacific C. glabellum Adams, 1868 from Japan. The evolution of a tubular shell marks the origin of Caecidae and likely correlates with a transition to an infaunal lifestyle (e.g., among corals and coral rubble or algae; Bandel 1996). However, interstitial habitats are very variable, differing with regards to the available space between the sand grains which influences the mobility, light intensity and therefore visibility and protection from predators. In the ‘Caecum glabrum-like’ microsnails, the morphological similarity among taxa (i.e., minute, slim shell, lack of ornamentation and coloration) likely correlates with a habitat shift into the mesopsammon and the consequent habitat restrictions of this special interstitial environment. ‘Regressive evolution’ leading to simplified and highly adapted body plans are typical for the mesopsammon (Swedmark 1968) and consequently the associated meiofauna is prone to cryptic speciation (Jörger et al. 2012; Meyer-Wachsmuth et al. 2014; Leasi et al. 2016).