Evaluation of the taxonomy of Helix cincta (Muller, 1774) and Helix nucula (Mousson, 1854); insights using mitochondrial DNA sequence data

Phylogenetic relationships of the Aegean Helix cincta and Helix nucula with congeneric species found in Greece were inferred using mitochondrial DNA sequence data. Twenty-three specimens from mainland Greece, Aegean Islands, Cyprus and North Africa were analysed, revealing that (1) H. nucula is monophyletic, (2) H. cincta from Greece and Cyprus is paraphyletic, and so questions arise regarding the taxonomy of this species, and (3) H. cf. cincta from Tunisia might be considered as a distinct evolutionary lineage. Moreover Helix valentini, an endemic species of Kalymnos Island group, is clustered within the lineage of H. c. anatolica, so supporting the synonymy of the two ‘species’, and the elevation of H. c. anatolica to species level. Hence, our results stress the need for a taxonomic reconsideration of H. cincta in the Aegean Sea, indicating that sequence data can prove useful in overcoming taxonomic issues at both species and subspecies level.


Introduction
The genus Helix is a large and diverse group of land snails that comprises more than 30 species in the Eastern Mediterranean region (Vardinoyannis 1994). They are among the largest European snails, which, apart from their scientific interest, are of increasing nutritional, commercial and economic importance, being widely consumed in Europe, particularly in Mediterranean countries. Helix cincta (Muller, 1774) and Helix nucula (Mousson, 1854) are two dark-lipped species that are present mainly in the Eastern Mediterranean region (Figure 1). Helix cincta is a relatively large, herbivorous land snail that can be found in a variety of different habitat types, such as rocky terrains (rocky grasslands, rocky shrublands and forests) and riparian vegetation (Schütt 1996). It is distributed through northern Italy, the Balkans, on the Aegean Islands and in western and southern Anatolia reaching Lebanon (Schütt 1996). More specifically, H. cincta has three subspecies; H. c. cincta (Muller, 1774) is found in the central Mediterranean (Dalmatic coast, and Italy), H. c. ambigua (Mousson, 1859) in Greece (continental Greece, Ionian Islands, and Northern Sporades), and H. c. anatolica (Kobelt, 1891) in the Eastern Mediterranean (East Aegean Islands, Turkey and Cyprus). Helix nucula, a medium-sized herbivorous species that inhabits sandy coasts, phrygana, semi-arid and arid areas, has a more restricted distribution. In the Aegean Sea, H. nucula is restricted to the island of Anafi (Cyclades) (Mylonas 1982) and southern Crete (Vardinoyannis 1994). It is also found in northern Africa (Zilch 1952;Kaltenbach 1959). The island of Crete is the only region where both H. cincta and H. nucula occur (not syntopically). It should be noted here that many land snails are considered to have been introduced to the island of Crete and the South Aegean islands from other Mediterranean regions, such as northern Africa (i.e. H. nucula, Cochlicella acuta, Theba pisana) (Welter-Schultes 1998).
The taxonomy of the two helicids is based on shell and genitalia features. In this work we did not consider records of the focal species reported in various public databases dealing with land snails of the region, namely Fauna Europaea, AnimalBase and the International Union for Conservation of Nature (IUCN), due to the large discrepancies among their approaches, coupled with the overall lack of documentation of their claims. The distribution of H. nucula mentioned in IUCN's Red List contradicts published data (Paget 1979;Maassen 1981;Frank 1997). In addition, there are two papers (Neubert 2000;Örstan et al. 2005) that provide erroneous distributional data for this species, as indicated by Schütt (2005), Triantis (2006) and Triantis et al. (2008). Finally, the digital image of AnimalBase documenting the presence of H. nucula in Turkey is not clear, and the authors of the website state that the specimen could actually belong to H. cincta.
Although morphological characters have been used to define the phylogenetic relationships and validate the taxonomy of terrestrial gastropods (Ponder and Lindberg 1997;Hausdorf 2002), their remarkable phenotypic diversity has led in several cases to controversial taxonomy (Parmakelis et al. 2003;Douris et al. 2007;Kornilios et al. 2009). In such cases, sequence data have repeatedly enabled researchers to formulate and evaluate alternatives to the morphological one-species 'hypothesis' and validate or refute the taxonomy based on morphology alone. In this study, using sequence data from two mitochondrial genes, we are exploring the phylogenetic relationships of H. cincta

DNA extraction, amplification and sequencing
All individuals included in this study are listed in Table 1. Locations of the sampling sites are shown in Figure 1. During sampling we tried to cover the whole distributional area of H. cincta and H. nucula in the Aegean Sea. For this reason we used all the samples available in the collections of the Natural History Museum of Crete (University of Crete). All specimens were assigned to species based on genitalia features and shell morphology. Genomic DNA was extracted from small pieces of foot, using the DNeasy Blood & Tissue extraction kit (Qiagen, Hilden, Germany). Doublestranded polymerase chain reaction (PCR) was used to amplify a partial sequence of the mitochondrial genes (mtDNA), encoding the large subunit ribosomal RNA (16S rRNA) and the cytochrome oxidase subunit I (COI). Primers 16SAr-l and 16SBr-h (Palumbi 1996) were used to amplify the 16S rRNA segment of approximately 530 base pairs (bp). The PCR cycle programme comprised an initial denaturation at 94 • C for 10 min, followed by 35 cycles of 1 min at 94 • C, 1 min at 50 • C, and 1 min at 72 • C, using single Taq DNA polymerase (Applied Biosystems, Foster City, CA, USA). The cycle ended with 10 min sequence extension at 72 • C. The universal primers LCO1490 and HCO2198 (Folmer et al. 1994;Hatzoglou et al. 1995) were used to amplify a COI segment of approximately 700 bp. The PCR cycle programme was as above, except the annealing temperature was set to 48.6 • C. Subsequently, a new set of specific internal primers were designed from sequences obtained with the universal primers. The specific primers Hgona_FCOI: 5 -TAATTGGHGGKTTTGGWAATTG-3 and Hgona_RCOI: 5 -GTRTTAAAATTTCGATCYG-3 were used to amplify a segment of approximately 380 bp. The PCR cycle programme was the same as previously described.
Products of the PCR were purified using the ammonium acetate DNA purification method. Single-stranded sequencing of the PCR product was performed using the Big-Dye Terminator (v3.1) Cycle Sequencing kit on an ABI377 automated sequencer, following the manufacturer's protocol. Primers used in cycle sequencing were the same as those used in PCR amplifications. In addition to H. nucula and H. cincta five Helix species (H. figulina, H. godetiana, H. lucorum, H. pomatia and H. valentini) were included in the analysis (see Table 1).

Alignment and phylogenetic analyses
The alignment of the concatenated 16S and COI sequences was performed with MAFFT (Katoh and Toh 2008) and visually corrected. Alignment gaps were inserted to resolve length differences between sequences, and positions that could not be unambiguously aligned were excluded. COI sequences were translated into amino acids before analysis, and stop codons were not found. Sequence divergences were estimated based on Tamura and Nei's (1993) model of evolution using MEGA (v.5; Tamura et al. 2011). Phylogenetic trees were constructed using neighbour joining, maximum parsimony, maximum likelihood and Bayesian inference methods. Neighbour joining analysis was performed using MEGA and the Tamura-Nei model of evolution (Tamura and Nei 1993). Bootstrapping with 1000 pseudo-replicates was used to examine the robustness of clades in the resulting tree (Felsenstein 1985). Maximum parsimony analysis was performed with PAUP * (v.4.0b10; Swofford 2002), applying heuristic searches using stepwise addition and performing tree bisection-reconnection branch swapping (Swofford et al. 1996). Confidence in the nodes was assessed by 1000 bootstrap replicates, with random addition of taxa. Heuristic maximum likelihood searches were performed in RAxML Blackbox (Stamatakis 2006) with 100 bootstrap replicates. The models used for the Bayesian inference analysis [General Time Reversible (GTR; Rodriguez et al. 1990); + gamma (G) for 16S rRNA and Hasegawa-Kishino-Yano (Hasegawa et al. 1985) + gamma (G) for COI] were selected using jModelTest (v.0.1.1; Posada 2008) and the Akaike's information criterion (Akaike 1974). Bayesian inference was performed in MrBayes (v.3.2.1; Ronquist et al. 2012), with four runs and four chains for each run for 10 7 generations, and trees saved to file every 100 generations. To confirm that the chains had achieved stationarity, we evaluated 'burn-in' by plotting log-likelihood scores and tree lengths against generation number using Tracer (v.1.5.0; Rambaut and Drummond 2008). The -lnL stabilized after approximately 10 6 generations and the first 10 4 trees were discarded as a conservative measure to avoid the possibility of including random, suboptimal trees. A majority rule consensus tree ('Bayesian' tree) was then produced from the posterior distribution of trees and the posterior probabilities were calculated as the percentage of samples recovering any particular clade, where probabilities ≥ 95% indicate significant support (Ronquist et al. 2012).

Results
For the phylogenetic analyses, a data set of 38 concatenated (COI and 16S rRNA) sequences was used containing 10 H. cincta, two H. cf. cincta (from Tunisia, indicating that they are morphologically similar to H. cincta), 11 H. nucula, two H. figulina, three H. godetiana, two H. lucorum, two H. pomatia and three H. valentini). A total of 975 bp were aligned, of which 337 (34.5%) were variable sites and 264 (27%) were parsimony informative. A COI fragment was not successfully amplified from H. godetiana or H. valentini. The analysis of sequences revealed 23 unique haplotypes. Sequence divergence ranged from 0 to 40.8%. Table 2 shows the genetic distances among the major lineages of each species, as well as within each major lineage. All phylogenetic analyses (neighbour joining, maximum parsimony, maximum likelihood and Bayesian inference) produced trees with similar (not identical) topologies, in which several clades are statistically well supported. Unweighted parsimony analysis produced more than 10,000 equally parsimonious trees with a length of 768 steps (consistency index, CI = 0.608; retention index, RI = 0.721). The large number of equally parsimonious solutions was due to terminal branch swapping. Identical topologies were recovered for each of the four runs with the full data set, and the 50% majority-rule consensus tree of the 9 × 10 4 trees remaining after burn-in is presented in Figure 2.
According to the phylogenetic tree (Figure 2), H. godetiana branches off first, whereas the focal species of this study form a monophyletic clade with substantial statistical support. However, the phylogenetic relationships within the cincta/nucula clade remain unresolved.
Helix nucula seems to be monophyletic, consisting of four distinct lineages with unresolved relationships. The first two lineages, with very good statistical support, comprise six specimens from south-central and west Crete, and two specimens from Libya, respectively. The other two lineages include two specimens from southeast Crete and one specimen from the island of Anafi.
Helix cincta appears to be polyphyletic, with the relationships of all the nominal H. cincta populations being unresolved. The H. cincta populations that seem to be quite distinct from the remaining ones are the H. c. anatolica originating from Cyprus and the east Aegean islands. Notably, the three specimens of H. valentini from three different east Aegean islands cluster inside the H. c. anatolica clade.
According to the phylogenetic tree of Figure 2, Helix cincta. H. cf. cincta, H. nucula and H. valentini form a single clade, pointing to the high phylogenetic affinity of these forms. Within this clade, H. cf. cincta from Tunisia branches off first. This can be considered as an indication that this lineage could be a different species. Unfortunately, no evidence can be provided to support or reject this scenario and the issue will remain unresolved until more samples and populations are analysed from North Africa.
In agreement with current morphological taxonomy, H. nucula seems to be a distinct monophyletic species. Within the H. nucula clade, two main lineages were statistically well distinguished. The first corresponds to the North African populations and the second to the populations from southwest Crete. Three other specimens (two from southeast Crete and one from Anafi Island) were also included in the H. nucula clade, but their phylogenetic position is unclear. It should be noted that the phylogenetic proximity of the lineages from the South Aegean area with lineages from North Africa (Libya) has also been depicted in the bat species Pipistrellus hanaki (Hulva et al. 2007) and the plant species Androcymbium rechingeri (del Hoyo and Pedrola-Monfort, 2008). This could be an indication that the South Aegean populations of H. nucula had their origins in the North African region, confirming the statement of Welter-Schultes (1998). However, it is impossible to support the claim that the species has been introduced in the islands of Crete and Anafi, until more specimens from North Africa are incorporated in the analyses.
The phylogenetic relationships among H. cincta lineages are more complicated. The only well-supported lineage of H. cincta includes specimens from the East Aegean islands and Cyprus (H. c. anatolica), as well as, three specimens of the local endemic H. valentini from Kalymnos island group (Telendos, Nera and Agia Kyriaki -East Aegean). The phylogenetic position of H. valentini in the tree (within the lineage of H. c. anatolica), the low genetic distance (3.2%) and the morphological similarity (shell morphology and reproductive system) of H. valentini with H. c. anatolica support the synonymy of the two 'species', with H. c. anatolica raised to species level as senior synonym. The other lineages of H. cincta are found in west Crete, Skyros, Kerkyra and Peloponnisos. The absence of a clear phylogenetic pattern indicates the inability of the present data to resolve the phylogenetic relationships among them. Consequently, although our results underline the need for more data to resolve their phylogenetic relationships and evaluate their taxonomic status, they also stress the need for a taxonomic re-evaluation within this group.
As a whole, the examination of mtDNA variation in the land snails of the species H. nucula and H. cincta may contribute to clarifying their taxonomic status. Phylogenetic information can now be added to the knowledge of their morphology and distribution, producing a more accurate taxonomy for these species.