Phylogenetics and species delimitations of the operculated land snail Cyclophorus volvulus (Gastropoda: Cyclophoridae) reveal cryptic diversity and new species in Thailand

Recent conceptual and practical advances in phylogenetic species delimitation have enabled progressively robust biodiversity studies. Delimiting species in widespread taxa is an intriguing problem; the edible operculated land snail Cyclophorus volvulus (Müller, 1774) is a good example since it shows a high degree of shell and color variation along with a widespread distribution throughout Thailand. Taxonomic boundaries for C. volvulus were examined and clarified using a combined morphological and phylogenetic approach, the latter of which was based on both nuclear and mitochondrial gene sequences. Moreover, three species delimitation analyses were applied: Poisson tree processes (PTP), automatic barcode gap discovery (ABGD), and generalized mixed Yule-coalescent (GMYC). All phylogenetic trees revealed that C. volvulus was polyphyletic and comprised of three clades that coincided with their geographic distribution. The three species delimitation analyses concurred with the phylogenies and formed at least three groups. According to the results, C. volvulus s.l., as currently recognized, consists of three distinct species in Thailand: C. volvulus s.s., C. occultus sp. nov., and C. borealis sp. nov., which are described herein. Moreover, several of these highly distinct C. volvulus evolutionarily significant units (ESU) are likely to require urgent conservation attention.


Results
Molecular and phylogenetic analyses. The 16S ribosomal RNA (rRNA) (462 bp), cytochrome oxidase I (COI) (660 bp), 28S rRNA (585 bp), and concatenated dataset of all three DNA regions (1,707 bp) were aligned for the 29 individual C. volvulus samples, which represented 15 populations from Thailand, along with Cyclotus spp. and Leptopoma vitreum as outgroups. The concatenated dataset had a 44.55% GC content. Heterogeneity in base composition between sequences is known to affect phylogenetic inference, so the variations in base pair composition among sequences for the concatenated datasets were tested using χ 2 analysis (χ 2 = 29.54, degrees of freedom = 120, P = 1.000). The dataset included 609 variable sites and 443 parsimony-informative sites. The partition homogeneity test, performed using PAUP 4.0b10 with 100 replicates, showed no significant incongruence between the gene fragments (P = 0.097), while the K2P-distance between the taxa (without outgroups) ranged from 0.000 to 0.120 (inter-and intra-specific K2P-distances were 0.079 and 0.012, respectively).
The phylogenetic relationships of C. volvulus samples are shown in Fig. 1. There was significant statistical support for a phylogenetic tree where certain clades coincided with certain morphological groups. However, C. volvulus was polyphyletic; it was divided into three well-supported main clades that represented three geographical regions, i.e., Northern (clade A), East-Central (clade B), and Western (clade C) Thailand. The sister taxon for clade A was C. subfloridus, while the sister taxa for clade B and C were C. labiosus and C. fulguratus, respectively. Clade A was strongly supported by 100% neighbor-joining (NJ) and maximum-likelihood (ML) bootstrap replicates and 1.00 posterior probability (PP) of Bayesian inference (BI). Clade B was supported by 94% NJ and 96% ML bootstrap replicates and 1.00 PP of BI, and clade C by 100% of NJ and ML bootstrap replicates Geometric morphometrics of shell morphology. Clades A-C ( Fig. 1) were provided as an a priori group for each canonical variate analysis (CVA), which supplied a graphical display of the shape differences between them from the relationship between the two obtained canonical variate (CV) variables (Fig. 3). The shape variability explained 66.64% of the variance in the first canonical axis (CV1), with an eigenvalue of 0.7039, while shell shape variability explained 12.78% of the variance in the second axis (CV2), with an eigenvalue of 0.8358. Overall, the shell shapes of the three clades showed some overlap. However, there were also significant differences among clades based on permutation tests for Mahalanobis and Procrustes distances (P < 0.0001 to 0.0018 for C. volvulus). Accurate shell classification for each clade, along with cross-validation of discriminant function analysis, was 78% for clade A, 67% for clade C, and 63% for clade B (Fig. 1), with a 69% overall rate of reliability for identifying each C. volvulus individual. The overall significance tests were based on pooled variances (i.e., the average variance across all groups). www.nature.com/scientificreports www.nature.com/scientificreports/ A wireframe modification (Fig. 3) defined the shape change from the consensus configuration along the two CVs. For CV1, individuals located in the positive portion of the axis had a wider shell compared to those in the negative portion (represented by shifts in landmarks 6 and 8 for C. volvulus). CV1 showed differences in aperture, where individuals in the positive portion of the axis had a narrower shell aperture compared to the negative portion. The CV2 of the C. volvulus complex showed differences in shell shape, prominently defined by landmarks 7 and 9, which represent the side of the last whorl twist, compared to the other landmarks as the score decreased.
In the present study, most of the examined C. volvulus samples showed a high degree of similarity in shell size and protoconch and jaw patterns (Fig. 4). However, we found slightly different characteristics among clades A-C, such as shell and umbilicus color patterns. However, the possible type specimen was placed in clade C; its shell was most similar to C. volvulus.

Discussion
In this study, we confirmed the existence of three highly divergent C. volvulus s.l. lineages in Thailand; these genetically differentiated and geographically isolated populations were confirmed using multi-loci species delimitation. These lineages (Clades A, B and C in Northern, East-Central, and Western Thailand, respectively; Fig. 1) were first recognized in a previous study 6 , and further verified here with the analysis of new samples as highly divergent populations (Figs 1 and 2). Therefore, the three gene loci that we analyzed were generally consistent with the previous study on the relationship of the genus Cyclophorus 6 . The individual gene trees showed some differences in topology and support for each lineage. However, individual gene trees are often different from the species tree due to retention and incomplete sorting of ancestral polymorphisms, especially in the case of recently diverged species or populations 6,18-20 . All three species delimitation analyses revealed the three C. volvulus lineages identified in the phylogenetic trees ( Fig. 1) and clearly reflected three-to-thirteen distinct evolutionary significant units (ESUs; Fig. 2). PTP species delimitation analysis revealed additional candidate species within C. volvulus relative to the ABGD and GMYC methods (number of candidate species: PTP = 13, ABGD = 3, and GMYC = 3). Intra-specific K2P divergences for mitochondrial DNA in the three main candidate species ranged from 0.011-0.041, while the other Cyclophorus species ranged from 0.000-0.021. These species delimitation analyses were responsible for the high average intra-specific K2P genetic divergence observed for all candidate species. This relatively high intra-specific genetic divergence suggests that there are at least three separate species (clades A-C). Additionally, the results highlight differences in evolutionary patterns for the three nucleotide regions during inter-specific divergence and following demographic changes in candidate species. Notably, similar genetic distances were inferred to be indicative of distinct species in studies of other land snails 7,21 .
Generally, cryptic diversity is common among land snails that are widely distributed and/or fragmented in distribution 12 . As time passes reduced gene flow and/or a complete lack of genetic exchange between populations would promote new species due to allopatric processes. Indeed, morphology comparison among the three clades with their sister taxa underscored this supposition. Specifically, clade A was similar to C. subfloridus in aperture shape but different in shell shape, size, and color pattern, clade B was similar to C. labiosus in color pattern but different in aperture shape, and clade C was similar to C. fulguratus in shell size but different in color pattern. C. volvulus clades A-C ( Fig. 1) and their sister taxa may have recently undergone allopatric fragmentation within www.nature.com/scientificreports www.nature.com/scientificreports/ the large geographic distribution between the different localities 7,22 . Alternatively, they may be under similar selective pressure or niches that cause them to undergo convergent evolution 23 . Intriguingly, our PTP analysis provided a different number of candidate species compared to ABGD and GMYC methods. This discrepancy may demonstrate that neutral genetic variation in the snail is influenced by the environment. Compared to the three consistent candidate species identified with all delimitation methods, the remaining ten PTP-derived candidate species exhibited relatively little difference in shell morphological characteristics (e.g., size, shape, and color pattern). We conclude that the intra-specific genetic divergence for these additional PTP-identified species is too low to currently consider them as distinct species 7,21 , but the relatively high genetic divergence indicates they may split into different species in the future.
Although the external C. volvulus s.l. morphology among the three clades was highly similar, there were slight differences. Shell shape geometric morphometrics revealed significant divergence among these groups as well as considerable overlap in the shell shape, results that support their cryptic morphology. The possible type specimen, considered most related to C. volvulus in shell shape, was placed in clade C. In addition to the high degree of external morphology similarity, the genital systems of these groups also could not be used for classification and identification 6,11 .
In this study, we provided a robust provision for the previously reported relationships and conservation management of Cyclophorus 6 . We also identified priorities for ongoing and future investigations. Further, an extensive taxonomic revision of this genus would be valuable. The evidence presented in this study strongly suggests that the  Radula and jaw. Taenioglossan radula, each row contains 7 teeth with formula 2-1-1-1-2. Central tooth symmetrical semicircular shape composed of 5 denticles; central cusp well developed, largest and pointed tip; inner lateral cusps small, triangular shape with pointed cusp; outer lateral cusps very small with dull cusps. Lateral and inner marginal teeth tricuspid consisted 3 cusps; central cusp large, convex shape, and flanked with smaller and pointed head of one inner and one outer lateral cusps. At the edge of ribbon, outer marginal teeth bicuspid consisted 2 cusps; inner cusp small with pointed head; outer cusp very large triangular shape with pointed (Fig. 4D). Description. Shell medium, thick, solid with low conical shape, apex acute; whorls 5 to 6, convex; suture deep and wide. Periostracum exhibits color variation from dark to pale brown. Protoconch with thin narrowly discernible transverse ridges (Fig. 4B). Last whorl rounded, enlarge and with whitish colour. Shell colour usually light to dark brownish colour, sometime dark zigzag streak present; on periphery with narrow dark spiral band; lower periphery with wide striated stripe surrounded umbilicus. Aperture circular, slightly oblique; lip little expanded and reflected with white colour. Umbilicus widely open and deep. Operculum corneous, multi-spiral and little concaved center (Fig. 5C).

Radula and jaw.
Central tooth big and semicircular with 5 denticles, large at center and smaller arrangement of the other. Lateral teeth bicuspid, endocone small with sharp cusp, ectocone large. Marginal teeth bicuspid, have different arrangement into two layers of inner marginal teeth locates the same level of central tooth, the large sharp inflated ectocone, the outer marginal teeth locates at the edge of the ribbon with a little lower position from lateral teeth, with almost the same shape as inner marginal teeth. General structure of jaw similar to that of C. volvulus s.s. (Fig. 4H) The difference characters are sculptured with thicken and strong longitudinal parallel ridges and connected to each other with very strong parallel slanting transverse ridges (Fig. 5C).
Distribution. This new species is known from several localities in the northeastern Thailand. The application of the taxon name in previously publications by Nantarat et al. 6 as the C. volvulus s.l. are here reconsidered as the new species.
Remark. This new species can be distinguished from C. volvulus s.s. and C. borealis sp. nov. by having uniform brown to brownish shell colour, with light brown zigzag streaks. Aperture lip is a little bit more oblique than C. volvulus s.s. Shell pattern around umbilicus with wide spiral band and light brown zigzag streaks almost umbilicus area. This new species also supported by phylogenetic tree (Fig. 1) and geometric morphometrics (Fig. 3). www.nature.com/scientificreports www.nature.com/scientificreports/ Measurements. Shell height: 24-26 ± 0.04 mm. Shell width, ranges 29-37 ± 0.07 mm.
Description. Shell medium, funnel-shaped, umbilicate, solid, smooth surface, brownish, apex acute; whorls 5, convex; suture deep and wide. Periostracum thick corneous with brownish colour. Protoconch with thin closely discernible transverse ridges (Fig. 4C). Remain whorls with only thin irregular growth lines. Last whorl rounded, expand, with white colour background. Shell colour brownish, with dark zigzag streak present; on periphery with narrow dark brown spiral band; on lower periphery with narrow striated stripe surrounded umbilicus. Aperture circular, slightly oblique; lip little expanded with white colour. Umbilicus widely open and deep. Operculum corneous, multi-spiral and little concaved center (Fig. 5D).

Radula and jaw.
Central tooth semicircular 5 denticles with large center and smaller arrangement of the other. Lateral teeth bicuspid, endocone small with sharp cusp, ectocone large. Marginal teeth bicuspid, have different arrangement into two layers of inner marginal teeth locates the same level of central tooth, the large sharp inflated ectocone, the outer marginal teeth locates at the edge of the ribbon with a little lower position from lateral teeth, with almost the same shape as inner marginal teeth (Fig. 5F).
General structure of jaw similar to that of C. volvulus s.s. The difference characters are sculptured with strong longitudinal parallel ridges, connected to each other with solid parallel curving transverse ridges (Fig. 5I).
Distribution. This new species is known from the northern Thailand. The application of the taxon name in previously publications by Nantarat et al. 6 .
Remark. This new species can be distinguished from C. volvulus s.s. and C. occultus sp. nov. by having uniform brown to dark brownish shell colour, with dark brown zigzag streaks. The last whorl is relatively larger than C. volvulus s.s. and C. occultus sp. nov. Shell height is relatively more than the other two species. Shell pattern around umbilicus with solid band outside and light brown spiral line with light brown zigzag streaks about ¾ of umbilicus. This new species also supported by phylogenetic tree (Fig. 1) and geometric morphometrics (Fig. 3).

Conclusion
Morphologically cryptic species among land snails are an important problem for taxonomists, ecologists, and biogeographers. Study of these species provides valuable information for communicating the unique properties of separate evolutionary lineages and conservation management. The presented phylogenetic interpretation of C. volvulus s.l., based on mtDNA (16S rRNA and COI) and nuDNA (28S rRNA) sequence data, was conducted using multiple approaches, including molecular phylogeny, species delimitation, geometric morphometrics, and morphological data. All methods supported the presence of at least three geographically distinct candidate species in C. volvulus s.l. over broad regions (clades A, B, and C from Northern, East-Central, and Western Thailand, respectively). We clarified the species boundaries of C. volvulus and the two newly described species. We unveiled cryptic diversity in a biodiversity hotspot: the identification of these (at least) three cryptic lineages suggests additional species diversity. The new species may comprise the portion of the potential niche that is actually occupied, given that historical and biotic factors pose further restriction on new species' distribution 24 . Differences among new species' niches could reflect their evolution, an actual change in the potential niche, and/or a mechanism that results in isolation and speciation. However, these differences could also reflect morphological plasticity, where some PTP candidate species possess the same potential niche despite ostensibly expressing a different niche. This study also provides useful information to evaluate ESUs for conservation, especially in some parts of Northeast and East-Central Thailand where the snail populations have declined dramatically. Understanding the high level of lineage diversity in an operculated land snail in Thailand, is within an area designated as one of the biodiversity hotspots and a world conservation priority 25 .

Methods
Taxon sampling and morphology. Cyclophorus volvulus s.l. was collected from 15 populations (Tables 3   and 4). Species determination was based on the publications of Kobelt 1 , Reeve 3 , Müller 10 , and Benthem Jutting 26,27 . The shells were subsequently compared with the relevant type specimens and reference collections, including those at the Natural History Museum of Denmark (Zoological Museum), University of Copenhagen, Denmark (ZMUC), The Natural History Museum, London (NHMUK), and Chulalongkorn University, Museum of Zoology, Bangkok, Thailand (CUMZ). Living specimens were preserved. The foot tissues were fixed and preserved in 95% (v/v) ethanol, while the remaining parts of each specimen were preserved in 70% (v/v) ethanol for morphological study. In total, 51 specimens of C. volvulus and related species were examined. Adults shell were measured with Vernier calipers to the nearest 0.1 mm and photographed with a digital camera.
Details of the shell surface, protoconch, jaw and radula morphology (form and formula of radula teeth) were observed using scanning electron microscopy (SEM; JSM-5410 LV) at the Scientific and Technological Research Equipment Centre (STREC), Chulalongkorn University. DNA extraction, amplification and sequencing. Genomic DNA was extracted from the foot tissue using a DNAeasy Tissue Kit (QIAGEN Inc.). The mitochondrial 16S ribosomal RNA (16S rRNA) region and cytochrome c oxidase subunit I (COI) gene region were PCR amplified using the 16sar and 16sbr 28 and the LCO1490 and HCO2198 29 primer pairs, respectively. The nuclear 28S ribosomal RNA (28S rRNA) region was PCR amplified using primers 28SF4 and 28SR5 30 . The thermal cycling conditions for amplification of the 16S rRNA and 28S rRNA were 2 min at 94 °C, followed by 36 cycles of 30 s at 94 °C, 30 s at 50 °C and 90 s at 72 °C, and then a final 5 min at 72 °C. That for the amplification of COI was performed in the same condition except that the annealing stage was changed to 2 min at 42 °C and the extension time to 2 min. Products were resolved and visually checked by 1% (w/v) agarose gel electrophoresis with SYBR Safe staining and blue light transillumination. A www.nature.com/scientificreports www.nature.com/scientificreports/ QIAquick purification Kit (QIAGEN Inc.) was used to purify the PCR products before being sent for commercial automatic cycle-sequencing at Macrogen, Inc. (Seoul, South Korea).
Phylogeny inference. The mitochondrial DNA (approximately 462 bp of 16S rRNA and 660 bp of COI) and nuclear DNA (approximately 585 bp of 28S rRNA) sequences were edited and aligned using MUSCLE version 3.6 31 , a sub-program within the MEGA 6.01 software 32 , and then improved manually. The unambiguous of 16S rRNA, COI (all codon positions) and 28S rRNA fragments were aligned across all samples. All base frequencies and molecular character statistics were calculated using MEGA 6.01 32,33 .   www.nature.com/scientificreports www.nature.com/scientificreports/ All sequences were checked for ambiguous nucleotide sites and saturation before being subjected to phylogenetic analysis. The saturation graphs were plotted using uncorrected pairwise distances for transition and transversion substitutions to envision saturation and identify the taxa responsible. In the cases of accuracy of the combined data suffered relative to the individual partitions, the concatenate datasets of mitochondrial (16S rRNA and COI) and nuclear DNA (28S rRNA) were checked using the partition homogeneity test in PAUP 4.0b10, using 100 replicates 34 . The alignment was partitioned by gene fragments. However, the analyses were run for each gene fragment and also for the concatenated dataset. jModeltest 2.1.7 35 was used to calculate and select a GTR + G model that was then used in all analyses.
Phylogenetic trees (Fig. 1) were constructed using neighbor joining (NJ), maximum likelihood (ML) and Bayesian inference (BI). The parameters, including the rate matrix and gamma shape parameter (α) of the gamma distribution (based on 16 rate categories) were estimated. The alignment was partitioned by gene fragment. The NJ analysis was performed using PAUP* v4.0b10 34 , while the ML analysis was performed via the RAxML Web servers 36 . Bootstrap resampling 37 with 1000 replicates was used to assess branch support. The BI analysis was performed using MrBayes version 3.2.5 38 , where the tree space was explored using four chains of a Markov chain Monte Carlo algorithm (MCMC) and the optimum substitution model found above. The Bayesian analysis was run for 10 million generations (heating parameter = 0.03), sampling every 100 generations and then discarded 25% of the trees with burn-in. Convergence was monitored by showing the average standard deviation of the split frequencies (between 2 runs) was below 0.1% 38 . Support for nodes was defined as posterior probabilities (PP).
Species delimitation methods. Three species delimitation methods (Fig. 2) were performed 39,40 , since the support across multiple species delimitation approaches is generally superior to a single method 40 . Thus, the species boundaries of C. volvulus were investigated based on the (i) Poisson Tree Processes (PTP) 16 , (ii) Automatic Barcode Gap Discovery (ABGD) 41 and (iii) Generalized Mixed Yule-Coalescent (GMYC) 42,43 methods. While ABGD detects 'barcode gaps' in the distribution of pairwise genetic distances, GMYC and PTP use an input tree as a model to delineate the PSHs in speciation and coalescent processes 40 . The ABGD approach was run using the web server (http://www.abi.snv.jussieu.fr/public/abgd/abgdweb.html) and the default parameters. The PTP and GMYC 43 methods were also performed on the concatenated 16S rRNA, COI and 28S rRNA genes.
The delimitated species tree was produced using the relaxed log-normal clock algorithm implemented in the BEAST v1.8.2 package 44 . The GTR + G model was applied to re-construct the tree for 1 × 10 7 generations with sampling every 100 steps. The MCMC output was examined by consideration of traces in Tracer 1.6 45 and analyzed with TreeAnnotator 1.7.4 using all trees. A PP limit of 0.5 with maximum clade credibility tree was set. The GMYC method was performed using both a single-and a multiple-threshold. The output tree was optimized using SPLITS v.1.0-19 package for R. The PTP method was implemented based on the output tree and run in Python using the Environment for Tree Exploration package 46 . Geometric morphometric analysis. Photographs of 108 shells, including the shell-only collections and the 42 individuals for which we assembled the DNA sequence data from, were randomly ordered in tpsUtil v.1.49 47 . The shell photos were then digitized in tpsDig v. 1.40 48 to capture the bi-dimensional coordination of 13 landmarks (Fig. 3A). Landmark coordinates for all specimens were evaluated in the software package MORPHOJ v1.06d 49 . All landmark data were superimposed by Procrustes Fit, with the data transformed into a common coordinate system to eliminate other variations, such as isometric size variation and orientation, except for shape from the data 47 . The allometric components/and the other effects of allometric shape variation were removed by using MORPHOJ v1.06d 49 . All statistical analyses were performed using MorphoJ 49 . Shell shapes were plotted against the centroid size which used as proxy of shell size to perform. Multivariate regression was used to discover any correlation among the centroid size and shape variable (p < 0.001; 10,000 replicates of permutation test) and its statistical significance. This allowed us to identify the modes of change that account for the highest degree of variation and provided linear numerical values for further analysis. Canonical variance analysis (CVA) was used to test the shell shape variation with candidate species recognized by molecular analyses as an a priori group. Mahalanobis and Procrustes distances between pairwise comparisons were calculated for significant differences by the permutation test (10,000 iterations). The correct classification percentage of each paired species was evaluated using the leave-one-out cross-validation, which was the implement of discriminant function analysis.
Nomenclatural acts. The electronic edition of this article fit with the requirements of the amended International Code of Zoological Nomenclature (ICZN), and hence the new names contained herein are available under that the Code rules. The published work and the nomenclatural acts it encloses have been registered in ZooBank (http://zoobank.org), the online registration system for the ICZN. The LSID for this publication is: urn:lsid:zoobank.org:pub:458492DC-E3F5-4293-AE27-FFCAB7287FC5. The electronic edition of this paper was published in a journal with an ISSN, and has been recorded and is available from PubMed Central.