Diversity of Land Snail Tribe Helicini (Gastropoda: Stylommatophora: Helicidae): Where Do We Stand after 20 Years of Sequencing Mitochondrial Markers?

Abstract

Sequences of mitochondrial genes revolutionized the understanding of animal diversity and continue to be an important tool in biodiversity research. In the tribe Helicini, a prominent group of the western Palaearctic land snail fauna, mitochondrial data accumulating since the 2000s helped to newly delimit genera, inform species-level taxonomy and reconstruct past range dynamics. We combined the published data with own unpublished sequences and provide a detailed overview of what they revealed about the diversity of the group. The delimitation of Helix is revised by placing Helix godetiana back in the genus and new synonymies are suggested within the genera Codringtonia and Helix. The spatial distribution of intraspecific mitochondrial lineages of several species is shown for the first time. Comparisons between species reveal considerable variation in distribution patterns of intraspecific lineages, from broad postglacial distributions to regions with a fine-scale pattern of allopatric lineage replacement. To provide a baseline for further research and information for anyone re-using the data, we thoroughly discuss the gaps in the current dataset, focusing on both taxonomic and geographic coverage. Thanks to the wealth of data already amassed and the relative ease with which they can be obtained, mitochondrial sequences remain an important source of information on intraspecific diversity over large areas and taxa.

1. Introduction

Intraspecific diversity is an important source of information about the mechanisms responsible for the current species distributions. The phylogeographic perspective reveals geographic structuring, informs about past distribution range extensions and population size changes and is able to distinguish between different scenarios responsible for accumulation of diversity in a given area and lack thereof in another (e.g., [1]). For example, genetic diversity allows differentiating between species evolved in situ and recent immigrants. Uncovering how biodiversity emerges and is maintained thus requires combining both inter- and intraspecific perspectives.

Unfortunately, datasets covering simultaneously intra- and interspecific diversity patterns within a taxon above the genus level remain rare and highly incomplete (e.g., [2,3,4,5]). It is difficult to obtain a broader picture of the diversity and phylogeny of, say, a family-level taxon across its whole distribution range, sampling all genera, species and the intraspecific diversity. To this day, collections of sequences of mitochondrial genes remain, along with microsatellites, the main source of information on the intraspecific variation on large taxonomic and spatial scales (e.g., [5,6,7]). Describing internal diversity of its constituent species in terms of both intrapopulation genetic variation and geographic structure is a mammoth task. Comprehensive and reliable information about diversity and its distribution is only obtained by a dense and geographically balanced sampling of populations. Data can be often repurposed from earlier studies, but utilization of published data is often difficult: they may be published in a summarized form [8] and they are often not appropriately georeferenced (e.g., [2,9]. Moreover, the sampling may be spatially and taxonomically biased depending on the goals of the original studies and filling the gaps in coverage then promises a diminished chance of discovery, which may lower the motivation to do such work.

We present here a dataset that covers both the species-level and intraspecific diversity of a Western Palaearctic tribe of land snails, Helicini. The tribe comprises the largest land snails in the region, several of which are very common and represent a prominent part of the local faunas. Thanks to several recent studies compiled here, this group is currently among the most thoroughly studied land snail taxa considering the distribution of intraspecific lineages. Its parent family, Helicidae, currently represents probably the best sampled land snail family, with more publicly available sequences than much more diverse Clausiliidae or Camaenidae, whose members have also been the target of numerous studies. Intraspecific diversity, however, is covered in a substantial part of the species only within the tribe Helicini, the intraspecific data for other helicid clades are less comprehensive or outright missing.

We collated georeferenced mitochondrial sequence data from published sources combined with rich own unpublished data, together spanning nearly 20 years of research on this group across over 80 species (according to https://www.molluscabase.org (accessed on 25 December 2021), but note that the taxonomy of this group is still not fully settled) and totalling to 2566 analysed individuals. We review what these data revealed about the taxonomy of Helicini and phylogenetic relationships between taxa, but the main focus is a qualitative comparison of the intraspecific diversity between species and, in particular, between different regions. For the first time we are able to compare the geographic patterns of distribution of intraspecific lineages across the whole group, including species broadly distributed as well as those with restricted ranges, by putting side-by-side species from different clades and regions. As we aim to provide a comprehensive primer on this model group for anyone interested in its diversity or in reusing the data further, we thoroughly discuss the strengths and weaknesses of the dataset, focusing on the gaps in taxonomic and geographic coverage.

2. Materials and Methods

2.1. The Model Group

The members of Helicini are medium- to large-sized land snails (greatest shell dimension ca 2–6 cm, e.g., [10,11]). The tribe is naturally distributed in the Western Palaearctic (ca. 20–54° N, 0–53° E, possibly extending further up to 68° E) and contains around 85 currently accepted species (see below). Neiber and Hausdorf [11] estimated the minimum crown age of the tribe to ca 31 Mya. This dating may be disputed because the placement of the fossil used for calibration is ambiguous (see [12]), but the estimates it yields are compatible with other lines of evidence [13,14] and the Oligocene age of the group is likely given the helicid fossil record [15].

The first sequence of a partial mitochondrial gene of a Helicini specimen (AF126144) was published in GenBank on 21 April 1999. It was a 372 bp fragment of 16S rRNA gene from Helix pomatia, with no locality given and misidentified as Helix lucorum, which was used as an outgroup in a phylogeographic analysis of another helicid species [16]. Other early uses of Helix mitochondrial sequences focused on identification of processed snail meat [17]. Manganelli et al. [18] were the first to provide phylogenetic insight into the systematics of Helicini using mitochondrial data, when they found indications that its type genus Helix, as usually delimited at that time, was polyphyletic. Helix pomatia was the subject of the first study focused on a Helicini taxon, presented at the World Congress of Malacology in 2007 [19], but the results were published only years later [20]. In 2012, Kotsakiozi et al. [21] published an analysis of the genus Codringtonia, which became the first published comprehensive molecular phylogenetic treatment of any Helicini genus. The foundations for the phylogenetics of the whole tribe were laid between 2015 and 2016, when its internal relationships as well as its position within the family Helicidae were explored [11,22,23,24]. Since then, further studies involving representatives of the tribe appeared. In 2019, Helix pomatia became the first Helicini species with a complete mitogenome sequence [25,26]. By now, the mitochondrial sequences still represent the bulk of existing genetic data for the Helicidae family (including Helicini). Nuclear sequence markers used to date consist mostly of the internal transcribed spacers 1 and 2 of the ribosomal rRNA cluster [20,27] and a gene for histone H3 (e.g., [24,28,29]). Genomic data are currently starting to be used [30] and the first draft genome of a helicid species has been published recently [31].

2.2. Data Acquisition

The dataset analysed here includes sequences publicly available from GenBank and our own as yet unpublished sequences accumulated since 2011. We also re-sequenced some DNA isolates analysed earlier (mainly in [23]) to obtain longer fragments. Only data available to us as of 31 December 2020 were included. The published data were collated from 33 peer-reviewed publications published between 2004 and 2021 [11,18,20,21,22,23,24,25,26,27,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54]. The Barcode of Life Data System (BOLD, https://www.boldsystems.org/ (accessed on 25 December 2021)) did not yield additional data, as only three Caucasotachea vindobonensis sequences available in BOLD were not represented in GenBank.

Four mitochondrial markers have been used so far for phylogenetic analyses involving representatives of Helicini. Most commonly these were the genes for 16S rRNA (16S hereafter; in two different lengths of the amplified region: ca. 400 or ca. 810 bp) and cytochrome c oxidase subunit I (cox1; 655 bp). Kotsakiozi et al. [21] also included a part of the cytochrome c oxidase subunit II (cox2; 505 bp) and partial sequences of the 12S rRNA (12S) were used in some studies attempting to reconstruct relationships between species or genera [24,35,37]. Here, we also successfully tested the amplification and sequencing of the 361 bp part of cytochrome b (cytb), used by [28] in a study of the helicid subfamily Ariantinae, on several samples across the diversity of Helicini. We also sequenced the 3′ half of the cox1 gene and the span between the cox1 and 16S genes, including the tRNA-Val gene, in representatives of major lineages within Helicini (as in [41]). Additional data extending beyond the five loci above come from transcriptome sequencing [25]. Five individuals (1 Caucasotachea, 4 Helix species) were analysed and the data also contained partial sequences of mitochondrial protein coding genes. We visualized the availability for different loci by plotting the distribution of the sequence data along the mitogenome of H. pomatia.

The new data were largely produced using the primers listed in

Table 1. List of primer pairs used for PCR amplification.

Forward Primer	Reverse Primer	Notes on Use	References
cox1
LCO1490: 5′-GGTCAACAAATCATAAAGATATTGG-3′	HC02198: 5′-TAAACTTCAGGGTGACCAAAAAATCA-3′	most samples	[55] *
p-cox1-f: 5′-TCGGGACGGGTCTCTCTTTG-3′	HC02198: 5′-TAAACTTCAGGGTGACCAAAAAATCA-3′	a few Helix pomatia samples	forward: this study
COI-vind-f: 5′-TACTGTTTGGTGTTTGATGTGG-3′	COI-vind-r: 5′-ACAACATAG TAATTGCCCCAGC-3′	Caucasotachea vindobonensis	[40]
LCO1490: 5′-GGTCAACAAATCATAAAGATATTGG-3′	H2198-Alb: 5′-TATACTTCAGGATGACCAAAAAATCA-3′	a few samples	reverse [32]
g-COX1-f: 5′-TGGGACAGGTTTATCGTTACTG-3′	H2198-Alb: 5′-TATACTTCAGGATGACCAAAAAATCA-3′	3′ part of Helix godetiana sequence	forward: this study
COI_OK1F: 5′-TTGTWACTGCYCAYGCRTTTG-3′	HC02198: 5′-TAAACTTCAGGGTGACCAAAAAATCA-3′	used for some samples that failed to amplify with the standard primer pair, mostly from museum ethanol material	forward [42]
LCO1490: 5′-GGTCAACAAATCATAAAGATATTGG-3′	COI_OK3R: 5′-AAAGGTGGRTAAACAGTYCANCC-3′		reverse [42]
16S rRNA
16Scs1: 5′-AAACATACCTTTTGCATAATGG-3′	16Scs2: 5′-AGAAACTGACCTGGCTTACG-3′	most samples	both [56] **
Scs1-p: 5′-GAATTACCTTTTGCATAATGGA-3′	Scs2-p: 5′-GAAACTGACCTGGCTTACG-3′	a couple of Helix pomatia samples	this study
16Scs1: 5′-AAACATACCTTTTGCATAATGG-3′	16S_MN3R: 5′-GCTACCTTTGCACAGTCAGWG-3′	the 16S sequence in two parts, mainly used for older or improperly preserved samples	reverse [57]
16S-F: 5′-CGGCCGCCTGTTTATCAAAAACAT-3′	16S-R: 5′-GGAGCTCCGGTTTGAACTCAGATC-3′		[58] ***
16S-F: 5′-CGGCCGCCTGTTTATCAAAAACAT-3′	16S-Helcentr-R: 5′-AAGYTTCTAGGGTCTTCTCGTCT-3′	the 3′ half of the 16S sequence in three parts, used for particularly fragmented templates	reverse [23]
16S-Helcentr-F: 5′-AGACGAGAAGACCCTAGAARCTT-3′	16S-R: 5′-GGAGCTCCGGTTTGAACTCAGATC-3′		forward [23]
16S-Helinter-F: 5′-GTACYYTGACTGTGCAAAGGT-3′	16S-Helinter-R: 5′-CTAGTCCAACATCGAGGTCAC-3′		[23]
16S-F: 5′-CGGCCGCCTGTTTATCAAAAACAT-3′	g-centr-R: 5′-AGACAGTTACCGCCCATGCT-3′	the 3′ half of the 16S sequence in three parts, for an old sample of Helix godetiana	reverse: this study
g-centr-F: 5′-AGCATGGGCGGTAACTGTCT-3′	16S-R: 5′-GGAGCTCCGGTTTGAACTCAGATC-3′		forward: this study
g-inter-F: 5′-TGGCCCATGATTGGGGTCTA-3′	16S-Helinter-R: 5′-CTAGTCCAACATCGAGGTCAC-3′		this study
cox2
COII-F: 5′-AAATAATGCTATTTCATGAYCAYGC-3′	COII-R: 5′-GCTCCGCAAATCTCTGARCAYTG-3′		[59]
12S rRNA
12SGast_fwd2: 5′-AGTGACGGGCGATTTGT-3′	12SGast_rev3: 5′-TAAGCTGTTGGGCTCATAAC-3′	most samples	[35]
12Sam: 5′-AACTAGGATTAGATACCCCAYTAT-3′	12bm: 5′-CGAGAGTGACGGGCGATTTGT-3′	a few samples, shorter amplicon than with the other pair	[60]
cytb
UCYTB151F: 5′-TGTGGRGCNACYGTWATYACTAA-3′	UCYTB270R: 5′-AANAGGAARTAYCAYTCNGGYTG-3′		[61]

* A shorter version of the reverse primer [60] may also be used. ** Combinations 16S-F+16Scs2 and 16Scs1+16S-R were also used in some cases. *** Alternatively, primers 16Sar+16Sbr [62] may also be used.

as the first for each locus. Other primer combinations were employed for amplification and sequencing when the standard combinations failed. For example, an incomplete sequence was originally obtained for the cox1 fragment in Helix godetiana, so the rest of the fragment has been amplified and sequenced with a forward primer specific to that sequence and H2198-Alb as the reverse primer.

PCR conditions varied depending on primer pair, the polymerase used and quality of the DNA isolate. Lately we used the following protocol as it appeared most efficient with the best PCR and sequencing results. We run PCRs in 20 μL volume containing 1 μL of the DNA isolate, 1.5 mM MgCl₂, 1X Platinum II PCR buffer (Invitrogen), 0.2 mM each dNTP, 0.2 μM each primer and 0.16 μL Platinum™ II Taq Hot-Start DNA polymerase (Invitrogen). The PCR cycle was set to 2 min at 94 °C and 35 cycles of 15 s at 94 °C, 15 s at 50 °C, 15 s at 68 °C.

Amplification and sequencing of the 16S locus was without problems with the 16Scs1 and Scs2 primer pair using the protocol presented here. The primers LCO1490 and HC02198 worked well for cox1 in most cases, except for Helix lutescens and some Helix schlaeflii lineages, where PCRs were mostly unsuccessful. For 12S the 12SGast_fwd2 and 12SGast_rev3 primer pair worked well, but for cox2 most samples produced reads where the same sequence was visible in the background shifted by one nucleotide, meaning that either the primers or the PCR protocol were not optimal. With cytb (and occasionally 12S), we encountered sequencing difficulties with some samples due to stretches of 8–11 thymine bases in a row.

The chromatograms of all new sequences were visually checked for reading errors. Sequences downloaded from GenBank had in some cases to be edited. Remnants of primer sequences or poorly read sequence ends were trimmed and in rare cases (when multiple lines of evidence suggested that data from the given study were carelessly curated) highly suspect substitutions (in conserved positions of 16S or non-synonymous in cox1) were given ambiguity codes. Although it is possible that they are actually accurate, some sequences were omitted altogether (e.g., KR705008, KF114835, JQ240036) due to suspected sequencing errors and/or long branches when compared to closely related samples.

An infrequent, but existing issue encountered with mitochondrial markers is the amplification of nuclear pseudogene sequences (NUMTs; e.g., [63]). We identified several instances of probable NUMT amplification in the Helicidae family as well as in Helicini. In Helix ceratina and one Cepaea nemoralis (Linnaeus, 1758) from Italy, we initially obtained two different sequences for cox1 by varying reagent concentrations and cycling conditions for the PCR with the same DNA isolate. We identified the genuine cox1 by amplifying and sequencing the region spanning from the standard cox1 fragment into 16S (H. ceratina [41]) or by comparison with 16S phylogeny (C. nemoralis). In another C. nemoralis individual from Italy, belonging to the clade E of Grindon & Davison [64], the sequencing resulted in a chromatogram with numerous double-peaks, suggesting that we co-amplified two distinct fragments simultaneously. Other possible NUMT examples are MF564162 from Helix melanostoma and MF564169 from Eobania P. Hesse, 1913 [29], which do not contain unexpected stop codons but, unlike the 16S from the same specimens, fall outside the correct clades when included in phylogenetic analysis.

All sequences, both newly obtained and retrieved from published studies, are listed with their metadata in Table S1. The geographic coordinates of sampling sites are given with varying precision, depending on how precisely the original location was known. In some cases, museum samples or published sequences have been used where the locality has been only verbally described, sometimes vaguely (for example, providing only the name of the closest settlement) or referring to places that we could not trace. In a few cases the published coordinates were corrected to correspond with the locality description. The museum lots indicated as vouchers include any shell material of the given species collected at the same site on the same occasion. Occasionally, these do not include shells of the sequenced individuals, typically when only small juveniles were found alive, which were preserved whole and directly used for DNA extraction.

2.3. Phylogenetic Analyses

2.3.1. Outgroup Selection

Currently, the best supported hypothesis on the relationships within the subfamily Helicinae assumes that the subfamily consists of two major clades, with centres of diversity in the western and eastern Mediterranean, respectively [11,25,27]. The tribe Helicini equals the eastern clade; the western clade consists of Allognathini Westerlund, 1903, Maculariini Neiber, Korábek, Glaubrecht and Hausdorf, 2021 and Thebini Wenz, 1923 [27] (but see [65] for an alternative system). Being sister to Helicini, the western clade is a natural outgroup choice.

For the outgroup, we collated only the 16S and cox1 data, because the other loci are available for only a few species. To avoid long branches as much as possible we employed a broad sampling of the outgroup by including most of the species recognized by the most recent revisions; more than one sample per species were included for Cepaea nemoralis (due to its very high intraspecific diversity) and for Cornu Born, 1778 (as the taxonomic splitting to species level is uneven across the genus). We mostly relied on published sequences [11,12,13,16,20,22,23,27,29,49,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82]. The respective GenBank accession numbers are listed in Table S2. Some species could not be included due to lack of available samples or sequences at the time of dataset collation (see [14] for additional data published meanwhile).

2.3.2. Alignment

The alignment of cox1, cox2 and cytb was straightforward since there were no indels when aligned with MAFFT 7.471, G-INS-i algorithm [83]. In contrast, alignment of the two rRNA genes is problematic due to the stem-and-loop structure of the rRNA, where the positions within loops may be non-homologous between distantly related sequences. Homology may be expected more safely among related species and these regions may also hold valuable phylogenetic information at the finest phylogenetic levels. In order to give priority to correct alignment of the loops among more closely related taxa, we performed the alignment in several steps. First, we aligned the sequences in each genus separately with MAFFT’s E-INS-i algorithm (default settings). The only exception was Cepaea, where we enforced the assumption of global homology by using G-INS-i. The sequences of the two Cepaea species are highly divergent and Cepaea hortensis (O. F. Müller, 1774) was underrepresented, which resulted in long unaligned blocks.

The resulting alignment blocks were checked and occasionally corrected when misplaced parts of incomplete sequences or apparently misaligned positions were observed at the end of aligned sequences. Then, we aligned these blocks within well-supported clades found in a recent analysis employing also the nuclear ITS2 (including partial 5.8S and 28S rRNA) data [27]: Otala+Loxana+Massylaea+Eobania+Gyrostomella, Cornu+Cantareus+Rossmaessleria, Allognathus+Hemicycla and Helix+Maltzanella. Within Thebini and Helicini, the resulting alignments were then aligned with each other and the remaining genera. In Allognathini there was an additional step of aligning all genera with the exclusion of Cepaea (see Figure 4 in [27]). We then aligned Thebini with Maculariini, these two with Allognathini and, finally, this complete outgroup with Helicini. Aligning the sub-alignments was done using the --merge option of MAFFT and the E-INS-i algorithm.

The data from the transcriptome sequencing were aligned by codons with MUSCLE [84] in MEGA 7 [85]. For phylogenetic analyses, only the cox1, cox2, cytb genes recovered from the transcriptome were used.

2.3.3. Maximum Likelihood Analysis of Backbone Phylogeny

We first examined the species-level backbone phylogeny of Helicini using two datasets: one consisting of the partial cox1 and 16S sequences and including the outgroup (“outgroup” dataset), the other without outgroup and aiming to maximize the length of the alignment (“maxloci” dataset) by including full cox1 and partial cytb, cox2, 16S and 12S where available. The nuclear ITS2 alignment used in analyses of multiple concatenated genes by Neiber et al. [27] was also analysed to show potential differences to the mitochondrial phylogeny.

For each dataset, the partition scheme and substitution models were selected with ModelFinder in IQ-TREE 1.6.12 [86,87] after initially partitioning the data into three codon positions for each protein-coding gene and separate partitions for each rRNA gene. Maximum likelihood analyses were then run with IQ-TREE and branch support was assessed with standard bootstrap (1000 pseudoreplicates) and SH-aLRT [88] (1000 replicates).

The inferred position of the genera Theba Risso, 1826 and Eremina L. Pfeiffer, 1855 within the outgroup, although without support, conflicted with expectations based on previous phylogenetic analyses (e.g., [27]) and morphology (Eremina appeared sister to Allognathini, Theba to Rossmaessleria P. Hesse, 1907). This appears to be driven at least in part by uneven nucleotide composition in different genera of the outgroup (see Results in Section 3). We also repeated the analysis with constraints on the outgroup topology to ensure that it had no effect on the inferred root position of Helicini. In this constrained analysis, we enforced the monophyly of Thebini with the exclusion of Macularia Albers, 1850 (according to [27]) and assigned Theba as sister to Eremina (due to shared preference for arid habitats, thick digitiform glands with reduced terminal branches and a small protoconch that is darkly coloured in some individuals).

2.3.4. Complete Phylogeny of Helicini

The strength of the presented dataset is in the coverage of intraspecific diversity of a number of the Helicini species. In order to describe and visualise the intraspecific lineage diversity in a manner allowing for comparison between species and regions, a unified objective approach to the delimitation of intraspecific groupings is needed. To this end, we constructed a complete time-tree of the samples and defined species-level clades and intraspecific clades by applying common clade age thresholds across the whole tree. The focus here were the phylogenetic relationships within species and the relative timing of diversification in different clades, not the relationships between more distantly related species and absolute dating of the tree.

Phylogenetic analysis of the complete dataset is challenging for several reasons. The number of samples is high, the alignment short and there is a great variation in the sequence length due to missing data [89]. Many of the samples yielded identical or nearly identical sequences and the loci available are insufficient to resolve deeper nodes within Helicini, especially relationships between related species and between genera [11,23,27,37,41], so it would be difficult if not impossible to obtain the complete phylogeny in one analysis. We therefore followed an approach inspired by that of Upham et al. [90]. We first constructed a backbone phylogeny using a dataset containing representatives from all major clades within Helicini and all samples of species unassigned to any of these, then we analysed each clade separately and, finally, we combined the resulting trees to create a complete phylogeny of all samples. The following datasets were analysed: Helicini backbone (cox1, 16S, 12S), Caucasotachea (cox1, 16S), Codringtonia (cox1, 16S, 12S, cox2), Levantina (cox1, 16S, 12S), Helix (Pelasga) (cox1, 16S, 12S), Helix Anatolian clade (cox1, 16S, 12S), Helix Mediterranean clade (cox1, 16S, 12S) and Helix European clade (cox1, 16S, 12S). We have run the single clade analyses without outgroups, because while each of the datasets (except for the backbone) comprises a well-supported clade, their closest relatives are in all cases uncertain [23,27] (see Results in Section 3 and Figure S1) and rooting of the clades could be biased if the outgroup is chosen arbitrarily.

For each dataset, partition scheme and substitution models were selected by IQ-TREE; we initially partitioned the data into three codon positions for each protein-coding gene and separate partition for each rRNA gene, but in all cases the model selection suggested four partitions: three codon positions and the rRNAs. We used BEAST 2.6.3 [91] for the phylogenetic analysis. We linked the tree and clock model between partitions. Bayesian Skyline was used as a flexible tree prior suitable for combination of inter- and intraspecific data [92] and a relaxed lognormal clock with mean rate of 0.02 substitutions per site per million years [13] was assumed for all analyses. In some cases, the selected substitution model has been downgraded from GTR to TN93 due to convergence problems of some of the model’s substitution rates in preliminary runs.

We run analyses of all sub-alignments in two replicates for 70 million generations sampling each 10,000th. 28% generations were discarded as burn-in after checking that parameter estimates converged within these and effective sample sizes over 200 were reached in the post-burn-in. This resulted in 10,000 trees for each sub-alignment. An exception was the European clade of Helix, the largest dataset (1292 tips) containing large amounts of very closely related samples mostly from H. pomatia. Analysis of this dataset was run in four replicates for 21 million generations.

The results for each single-clade dataset were summarized on a maximum clade credibility (MCC) tree with mean node heights (trees available from the Dryad repository). There were large differences in support for the position of the root between the single clade analyses and the backbone BEAST tree, with resolved root position only in the latter analysis (except for Levantina and Codringtonia supported in both). In the case of the European clade of Helix, there was a strong support for a root between Helix lutescens and the remainder in the backbone analysis (Figure S3). Helix lutescens is the only species in the clade which does not live in the Balkans [10] and lives syntopically with other members of the clade without any sign of past or present hybridization (own observations). Therefore, we consider this rooting very likely; it also appeared in earlier analysis with an outgroup [37]. In addition, the root for Caucasotachea was fully supported in the backbone tree. Its position corresponded to that uncovered in an earlier analysis [24] and is likely based on conchological and geographic grounds. Finally, we considered more likely that the root of the Mediterranean clade is between H. ceratina (or H. ceratina+H. ligata complex) and the remainder of the clade than among the conchologically similar species with brown shell apertures [41], as the single-clade analysis suggested (although without statistical support). We therefore accepted the root positions as inferred in the backbone analysis, despite uncertainty in case of Pelasga and the fact that the high supports for the root positions from the BEAST backbone analysis (Figure S3) were not mirrored in the ML analyses. We excluded from the posterior of the backbone and single-clade analyses trees not conforming the respective root positions. MCC trees were calculated for these filtered posterior samples and used for constructing the complete phylogeny. We found no appreciable effect on the support values for species-level and intraspecific clades (which were the focus of our analysis, see below); the differences only concerned the basal relationships between species which were not statistically supported in either analysis. The filtering led to varying reductions in the number of the posterior trees: for example, while no trees were excluded for Codringtonia and Levantina, almost a half did not correspond to the assumed root position in Caucasotachea and only ca. 10% of the posterior trees conformed to the selected rooting with H. ceratina in case of the Mediterranean clade.

The MCC trees from single-clade analyses were then grafted onto the MCC tree from the filtered posterior of the backbone analysis in place of the respective clades. The tree heights were adjusted to the height of the most common recent ancestor (MRCA) of the respective clades in the backbone tree. Note that there were substantial differences in the heights between single clade analyses and the backbone analysis despite using the same average clock rate, which corresponds to variability in the inferred clock rates for major branches in the backbone tree.

2.4. Distribution Maps of Intraspecific Lineages

We recognized two levels of clades for the visualization of the distributions of mitochondrial diversity (“species-level” and “intraspecific”) to obtain comparable units for plotting the distribution. We set these common age thresholds in the complete phylogeny and recognized the (sub)clades, whose crown ages were younger and stem ages older than these thresholds. Note that due to variation in molecular clock rate within the tree and the assembly of the complete tree from several subtrees, these thresholds may in fact, correspond to somewhat different absolute ages in different clades and subclades.

The “species-level” threshold has been set only for plotting and does not reflect any taxonomic opinion. No such threshold, however, would fit all situations. For our purpose, the “species-level” threshold was set so that well-established species were not split into multiple “species-level” clades, but in some cases the resulting clades contain more than one currently recognized species (e.g., Helix lucorum and Helix nicaeensis; Levantina). The MRCA of Codringtonia parnassia was used to set the “species-level” threshold. The threshold defining the intraspecific clades delimits subclades corresponding to or finer than divisions used in earlier phylogeographic studies [38,40]. The threshold was set at the base of Helix pomatia clade F sensu Korábek et al. [38] and adjusted so that several lineages represented by 1–2 individuals became included in larger intraspecific groups.

The geographic distributions of intraspecific lineages of all “species-level” clades were mapped using the same map scale and projection for direct comparability.

2.5. Non-Native Populations

Several large helicids (C. nemoralis, C. aspersum (O. F. Müller, 1774), T. pisana (O. F. Müller, 1774), E. vermiculata (O. F. Müller, 1774); for Helicini, see [23,38,39,41,42]) have been subject to countless intentional or unintentional introductions beyond their natural range limits. For meaningful biogeographic considerations, samples originating from such non-native populations must be identified and filtered out. We considered a sample to be from a native population when it was taken in a part of the species’ range where there is no a priori reason to doubt its natural occurrence. Populations from parts of the range which are known or suspected by us to be a result of introductions were labelled as non-native. The status of samples from near the tentative limits of natural distribution of H. pomatia and H. lucorum was given as unknown.

3. Results and Discussion

3.1. Mitogenome Representation

The mitogenomes of Helicidae sequenced so far are from 14,050 (Cornu aspersum) to 14,795 (Theba pisana) bp long [73,93,94]; the three sequenced mitogenomes of H. pomatia were very similar and 14,070–14,072 bp long [25,26]. For most samples, only about 10% of the mitogenome length or less was sequenced (Figure 1). In the presented dataset, only 33 out of the total of 2566 analysed individuals were sequenced for all four focal genes (cox1, cox2, 16S and 12S). 1421 were analysed for cox1 (max. fragment length 1506 bp, min. 133 bp, median 655 bp), 122 for cox2 (max. 658 bp, min. 360 bp, median 505 bp), 2313 for 16S (max. 851 bp, min. 122 bp, median 776 bp) and 121 for 12S (max. 675 bp, min. 240 bp, median 640 bp). The cox1 fragment is generally under-represented in comparison to 16S, but in Caucasotachaea many individuals were sequenced for cox1 only.

3.2. Dataset Coverage

The coverage of the distribution range of Helicini by currently available mitochondrial sequence data is highly uneven (Figure 2). The most densely sampled region are Central Europe and the Balkans. Towards the east, the sampling gradually becomes sparser. Those less covered areas include vast expanses of Eastern Europe with only few broadly distributed lineages of Helicini, as well as Anatolia and the Middle East where many endemic lineages occur.

The phylogenetic coverage is biased as well. Some species (especially H. pomatia, C. vindobonensis) are very well sampled across most of their range, but there are also species represented by samples from a single locality and genera like Levantina and Isaurica are only incompletely sampled at the species level. In fact, H. pomatia alone makes for 40% of the sequenced individuals and C. vindobonensis for additional 14%. This is partly due to studies dedicated to the phylogeography of these two species [38,40,51], but also due to their exceptionally large distribution range and high abundances in well accessible parts of Europe. In contrast, Helix salomonica from eastern Turkey, northern Iraq and Iran is an example of a broadly distributed species with very sparse coverage (three localities only).

In total, only 6% of the analysed individuals came from non-native populations and in further 4% the origin was classified unknown. However, in the extreme case of Helix cincta s. str., only one out of 30 analysed samples originated from its presumed native range.

3.3. Phylogeny

The phylogenetic analyses yielded, as expected, results similar to earlier studies. The available data do not provide a resolved mitochondrial phylogeny of Helicini (even with the “maxloci” dataset, Figure 3). The tribe Helicini is consistently recovered as monophyletic with both mitochondrial (Figure S1) and ITS2 data (Figure S2; see [27] for results of a concatenated analysis), but the only grouping between genera, supported unambiguously by the mitochondrial data, is the sister relationship between Helix and Maltzanella. The root position suggested by the analysis with outgroup (Figure S1; regardless of the outgroup topology) as well as the molecular clock analysis (Figure S3) is between Helix+Maltzanella and the rest of the tribe, but without unambiguous support. The supported groupings of species correspond to genera and the clades within Helix recognized by Korábek et al. [23]. Relationships among species are generally unresolved, but well-supported clades are found at the level of species and within species.

3.4. Diversity of Helicini

The tribe Helicini is unambiguously supported by molecular phylogeny, geographically well-defined and its delimitation is now agreed upon. The phylogenetic relationships among the genera within the tribe, however, are not resolved with the available mitochondrial data or by ITS2 data (Figure S2). In addition, the limited information content of the sequenced genes and missing data in the matrix, the evolutionary history of the group also contributes to the issue. Two genera are monotypic, in 3–5 others the crown group is young relative to the age of their branch. The situation is similar within Helix, as relationships among the main clades remain obscure. Furthermore, at least two of these clades (Mediterranean, European) seem to have undergone a rapid initial diversification, another factor contributing to poor resolution of the phylogeny.

3.4.1. Caucasotachea C. Boettger, 1909

In their morphology-based revision, Neubert and Bank [95] restricted the genus Caucasotachea to four species distributed in the Caucasus region, Alborz mountain range in Iran and along the southeast of the Black Sea by excluding species of the former subgenus Lindholmia P. Hesse, 1919. Molecular phylogenetic studies supported this split, but nevertheless led to changes in both the genus delimitation and species-level taxonomy. The very broadly distributed east-European Caucasotachea vindobonensis (C. Pfeiffer, 1828) has been transferred to this genus from Cepaea Held, 1838 [11,23,24,35] and the original four species were reduced to only two by Neiber et al. [47] following a detailed molecular genetic study.

The phylogeography of two species, C. vindobonensis (Figure S4) and Caucasotachea atrolabiata (Krynicki, 1833) (Figure S5), has been studied comprehensively [40,47,51]. The distribution of the former is well covered by sampling except for its eastern part from eastern Ukraine and south-western Russia to the river Volga and the Caucasus. In addition, its eastern range limits are not yet clarified. Due to the long distance between glacial refugia in the Balkans [40] and the eastern range extremes and the ubiquity of the species, C. vindobonensis may be a good model for studying how the genetic diversity becomes depleted with increasing distance from the postglacial expansion source in the absence of major dispersal barriers. For such a purpose, additional data from its eastern populations would be beneficial.

In C. atrolabiata, our new data from Abkhazia reveal that all major mtDNA clades are present, possibly in a parapatric pattern, in the western end of the Greater Caucasus, suggesting that its diversification centre lies here. To confirm the pattern, more samples from Abkhazia would be desirable to fill a gap in sampling.

Only a few sequences were available for Caucasotachea leucoranea (Mousson, 1863), which is distributed around the south of the Caspian Sea from Azerbaijan to Golestan in Iran (Figure S6).

3.4.2. Neocrassa Subai, 2005

Neocrassa has been distinguished first as a subgenus of Codringtonia [96], but the mitochondrial phylogeny indicates that it has to be separated at the genus level [23,27]. Its only species, Neocrassa neocrassa (Zilch, 1952), has a very limited range at the Greek–Albanian frontier (Figure S7). The available data revealed only shallow intraspecific divergences; no samples were available from Albania.

3.4.3. Isaurica Kobelt, 1901

The genus Isaurica is distributed in a small area in south-western Anatolia. The last revision [97] distinguished six species. However, Nordsieck [98] suggested that Isaurica callirhoe (Rolle, 1894) belongs, based on shell microsculpture, to Levantina and Amanica praecellens (Nägele, 1901) has been excluded from the genus by molecular phylogeny even earlier [23]. Of the remaining four species, Isaurica riedeli Subai, 1994 from north of Manavgat and near Akseki in the Antalya Province and Isaurica pamphylica Subai, 1994 from around Sütçüler and the Köprülü Canyon in the Isparta and Antalya Provinces of Turkey remain unsampled for molecular analyses. The phylogeny revealed samples of Isaurica schuetti Subai, 1994 as paraphyletic relative to Isaurica lycia (Martens, 1889), but the significance of this is unclear as only two individuals were analysed for each of these taxa (Figure S8). Introgression from I. lycia seems a plausible explanation.

3.4.4. Amanica Nordsieck, 2017

Amanica is a narrowly distributed monotypic genus separated from Isaurica after sequence data became available [23,98]. It only comprises Amanica praecellens (Nägele, 1901) from the Hatay province of Turkey (Figure S9).

3.4.5. Levantina Kobelt, 1871

Levantina comprises more than 20 currently accepted rock-dwelling species broadly distributed in the Middle East from Cyprus, Central Taurus mountains and the Levant to western Iran, Iraqi Kurdistan and south-western Arabia. The centre and origin of the present diversity of Levantina lies in eastern Turkey in the area south and south-west of Lake Van towards the Syrian Desert [42]. The species-level taxonomy is not resolved due to scarcity of samples from some regions and taxa on the one hand and variability of shell characters on the other hand. A thorough revision would require extensive new sampling. Intraspecific variability and its distribution are virtually unknown within Levantina. The genus has been recently newly delimited by excluding one former subgenus as a completely unrelated lineage and merging Assyriella P. Hesse, 1909 and Levantina (Laevihelix) Neubert, 1998 with the nominotypic subgenus [27]. In a follow-up study, the first mitochondrial phylogeny of the genus was presented [42], which forms the basis of the data shown here.

The western limits of the natural distribution of Levantina are uncertain. The populations of Levantina spiriplana spiriplana (Olivier, 1801) and L. spiriplana caesareana (Mousson, 1854) in the south-eastern Aegean and on Cyprus are most likely introduced [42], but there are two additional taxa from the western end of Levantina’s range, which have not been studied yet by molecular methods. One is Levantina rechingeri Fuchs & Käufel, 1936, known only from a few empty shells found on the slopes of the Kali Limni mountain on Karpathos island in the south-eastern Aegean [99,100,101]. The other is the above-mentioned Isaurica callirhoe (Rolle, 1894), which Nordsieck [98] reassigned to Levantina. It is known only from shells of the type series collected at an unknown location on the northern slopes of Akdağlar between Fethiye and Elmalı [102] in the very south-east of Anatolia.

Similarly, the range extent of Levantina in the east is not well documented. The known distribution extends roughly to Tehran [102], but Levantina longinqua (Schütt & Subai, 1996) has been described based on shells allegedly from “Hasrat Sultan Gebirge” south-east of Samarqand, Uzbekistan [103]. It remains known only from its type series and the type locality is not given precisely: Khazret-Sultan (Hazrati Sulton) at the Tajikistan border is the highest peak in Uzbekistan, but a broad area around it has to be considered. Finally, the true extent of distribution of Levantina in the west of the Arabian Peninsula is unclear and Levantina semitecta Neubert, 1998 from an unknown locality (probably in an area roughly between Jabal al-Lawz and al-Wajh in the north-west of Saudi Arabia) is known only from its two types collected in the 19th century [104].

In addition, the four problematic taxa above, no data are currently available also for Levantina mahanica Kobelt, 1910, described from near Lake Urmia and distributed south of it into Iraqi Kurdistan, and Levantina ninivita (Galland, 1885), described from near Mosul and recorded also near Cizre in Turkey [102]. In addition, samples from type localities would be desirable for Levantina guttata (Olivier, 1804) (Turkey, Şanlıurfa Castle hill) and Levantina thospitis (Schütt and Subai, 1996) (Turkey, between Bitlis and Baykan, Kermate/Alaniçi SW of Şetek/Ortakapı). They are likely conspecific with Levantina vanensis (Schütt and Subai, 1996) and Levantina mardinensis Kobelt, 1900, respectively, but because the intraspecific variation and significance of conchological characters in Levantina are poorly understood, the two pairs were not formally synonymized based on the available samples [42].

The majority of Levantina taxa, including the type species L. spiriplana, group in a large, broadly distributed clade with very short and unresolved branches at its base (Figures S10 and S11). This group is similarly aged as the Peloponnese radiation of Codringtonia (Figure S3; see also Figure 4 in [27]), but its range extends from southern Israel, Cyprus and the Cilician Taurus to the very south-east of Turkey. We presume that L. mahanica also belongs here as a potential close relative of Levantina kurdistana (L. Pfeiffer, 1862), which would extend the distribution of this clade up to Iraqi Kurdistan and western Iran.

A clade comprising Levantina djulfensis (Dubois de Montpéreux, 1840) and its relatives Levantina ghilanica (Mousson, 1876) and Levantina mazenderanensis (Kobelt, 1883) is distributed in the north-eastern part of the range of the genus (Figure S12). Levantina djulfensis has been reported also from an isolated area south of Siirt (Schütt and Subai 1996), a record whose identification may be worth a revision using molecular data due to its position in a region of high diversity of Levantina. Levantina ceratomma (L. Pfeiffer, 1856) (Figure S13) appears to be an isolated and well recognizable species, its samples from Iran west of Lake Urmia are nevertheless missing. Levantina escheriana (Bourguignat, 1864) is relatively broadly distributed (Figure S14). It is likely closely related to or conspecific with L. ninivita, judging from conchological similarity. It was found to be a sister clade of the three Arabian taxa Levantina asira Neubert, 1998, Levantina symensi Neubert, 1998 and Levantina asagittata Neubert, 1998, which are very closely related and for which samples were available only from their type localities (Figure S15).

The basal-most mitochondrial lineage in Levantina was recovered from a population sampled near the south-western end of Lake Van (Figure S16). We could not assign this sample reliably to any of the species accepted by Schütt and Subai [102]. The lineage is deeply divergent from the rest of the genus; however, the nuclear ITS2 data do not confirm the deep divergence between this sample and the remaining Levantina sequences (Figure S2).

3.4.6. Codringtonia Kobelt, 1898

Codringtonia comprises large rock-dwelling snails from Central Greece and the Peloponnese. The last taxonomic revision [96] distinguished, based on morphology, seven parapatric species. They are, except for Codringtonia parnassia Roth, 1855, mutually similar and very closely related. There is substantial geographic variation in details of shell shape and colour also within some of the species, suggesting geography as the main driver of diversification of Codringtonia. In the first molecular phylogeny, Kotsakiozi et al. [21] found support for the proposed classification with the exception of Codringtonia gittenbergeri Subai, 2005.

Our new data, which also include original material from the Subai’s revision [96], suggest a more complicated situation. Two samples from within the range of C. gittenbergeri as indicated by Subai [96] were identified as Codringtonia codringtonii (Gray, 1834) by Kotsakiozi et al. [21]. After excluding these, they concluded that C. gittenbergeri shares the same mtDNA clade with Codringtonia elisabethae Subai, 2005. However, the typical C. gittenbergeri shells including the holotype have an appearance similar to C. codringtonii, they are just darker, and the two paratypes we analysed also had mtDNA corresponding to C. codringtonii. Therefore, we consider C. gittenbergeri a junior synonym of C. codringtonii. Samples identified as C. gittenbergeri in Kotsakiozi et al. [21] were collected where Subai [96] reported C. elisabethae and likely indeed belonged to that species as the mtDNA suggests. The ranges of C. elisabethae and C. gittenbergeri, as indicated by Subai [96], adjoin and there is some overlap in the distribution of the corresponding mtDNA lineages.

After revising the status of C. gittenbergeri, C. codringtonii has apparently a disjunct range in the Peloponnese (Figure S17), formed by two areas separated by the range of Codringtonia eucineta (Bourguignat, 1857) (Figure S18). The sample from the northern Peloponnese carrying a haplotype of C. codringtonii was originally identified as C. eucineta by Subai [96] but that author already noted that there is a similarity to C. codringtonii with regard to shell characters.

The most diverse species within Codringtonia is clearly C. eucineta (Figure S18), which is also the most broadly distributed and most conchologically variable one. Samples of Codringtonia intusplicata (L. Pfeiffer, 1851) belong to a clade distributed eastward of C. eucineta (Figure S19) and two shallowly differentiated sister clades corresponding to Codringtonia helenae Subai, 2005 and C. elisabethae (Figure S20) occur even more to the east.

The clade uniting C. codringtonii, C. eucineta, C. intusplicata, C. helenae and C. elisabethae falls below the “species” threshold, so the divergences between these taxa are comparable to those often seen within other species of Helicinae. In the central Peloponnese, we have found multiple cases of discrepancy between the identification based on shell characters and the mtDNA lineage of the respective individual. Most involve C. helenae, where haplotypes of the clades characteristic for C. intusplicata and C. eucineta were found. In three cases these occurred in the same population together with haplotypes of the C. helenae clade. In addition, one individual identified as C. eucineta was found to have mtDNA belonging to a clade characteristic for C. helenae. A lineage from the C. intusplicata clade was found also in C. eucineta at one site in southern Aetolia. These discrepancies are mostly attributable to introgression as they occur at the contact between specie ranges, but incomplete lineage sorting seems possible for two lineages at the base of the C. intusplicata clade. Despite overlap in the ranges of the mitochondrial clades (Figure 4), we are not aware of syntopic occurrence of two currently recognized Codringtonia species except for a shared locality of C. helenae and C. intusplicata reported by Kotsakiozi et al. [21]. In light of our results, we doubt that this was indeed a case of coexistence of two separate populations.

There is a considerable phylogenetic diversity within C. parnassia (Figure S21), comparable to that within the clade uniting all other species of the genus. We have uncovered additional divergent mitochondrial lineages within this taxon on top of those reported by Kotsakiozi et al. [21] and slightly extended the known distribution to the north-east compared to Subai [96]. Various populations of C. parnassia also differ substantially in shell size and shape. Apparently, the taxonomic treatment is not comparable between C. parnassia and the rest of the genus and further revisions within the genus would be warranted.

3.4.7. Lindholmia P. Hesse, 1919

This genus comprises two recognised species [94], both of which were sampled and analysed. Lindholmia christophi (O. Boettger, 1881) is known from a small area in north-western Turkey near Artvin in the vicinity of Borçka and Ardanuç (Figure S22; the samples come from the latter locality). Lindholmia nordmanni (Mousson, 1854) has a distribution extending considerably more to the west than is the extent of the sampled sites (Figure S23), up to the west of the Yozgat Province of Turkey.

3.4.8. Maltzanella P. Hesse, 1917

In his revision, Schütt [105] recognized two species within this genus: Maltzanella dickhauti (Kobelt, 1903) and Maltzanella maltzani (Kobelt, 1883). We have found two lineages above the “species” threshold in M. dickhauti (Figure S24) from south-western Anatolia, suggesting that there may be additional diversity of mtDNA lineages yet to be uncovered. Maltzanella maltzani is known only from a small area near İzmir in western Turkey, but Korábek et al. [23] reported a single Maltzanella individual collected in the European part of Turkey in Kuru Dağı which yielded a mtDNA haplotype close to the M. maltzani sample (Figure S25). The specimen (SMF 342502) was conchologically more similar to M. dickhauti; apparently, the conchological diversity of the genus and its relationship to phylogeny and taxonomy is not yet sufficiently known.

3.4.9. Helix Linnaeus, 1758

Helix was the first helicid genus, in which molecular phylogenetics demonstrated the necessity of changes in its taxonomic delimitation [18,106]. Some of those changes were proposed even earlier on the grounds of genital system anatomy [107,108], namely the exclusion of Cornu aspersum and Cantareus apertus (Born 1778) from Helix. Later works included within Helix the following two genera: Tacheopsis C. R. Boettger, 1909 from north-western Anatolia [11,23] and Tyrrhenaria P. Hesse, 1918, endemic to Corsica [20]. The first nearly complete molecular phylogeny of the genus distinguished four major clades and four unassigned species within the genus [23]. Neubert [10] recognized only two subgenera based on morphology of the genital system, Helix and Pelasga P. Hesse, 1908, but whether Helix is monophyletic in respect to Pelasga remains still unclear due to unresolved relationships between major groups within the genus.

We resolve here the last remaining issue regarding what taxa should be included in Helix. Neubert [10] proposed to transfer Helix godetiana Kobelt, 1878 from islands in the southern Aegean (Figure S26) to Maltzanella. There are substantial conchological similarities between those taxa and only a very short fragment of 16S was available at the time, which did not refute that hypothesis [23]. We were now able to obtain the complete 16S and cox1 fragments from a dry museum specimen collected in the late 19th century (a syntype of Helix dacoronae Letourneux, 1884). The results (Figure 3 and Figures S1 and S3) show H. godetiana with full support as a member of Helix. Its precise position is unresolved, but the results suggest it could be either the basal-most Helix species or a sister clade to the subgenus Pelasga. Further analysis of this rare species is warranted.

Subgenus Pelasga P. Hesse, 1908

Species of the subgenus Pelasga are distributed from Greece and North Macedonia along the Mediterranean coast to Israel and Jordan, with one species extending eastwards to Iran. The diagnostic character of the group is the epiphallus, at least twice as long as the penis. The shells are, except for some large forms from the Levant, also very similar between species and their shape and sculpture is characteristic for the subgenus. Neubert [10] recognized six species, but this will probably need a future revision.

Of all major Helix clades, Pelasga remains the least sampled because the snails are buried in the soil when inactive, so live individuals can usually be found only for a limited part of the year [109] and/or shortly after rain. Therefore, even the type species of the subgenus, Helix figulina Rossmässler, 1839, is not well represented in our dataset (Figure S27), although it is very abundant in a large part of Greece. Similarly, Helix salomonica Naegele, 1899, a species with the eastern-most distribution, has a vast range in south-eastern Turkey (westwards at least to Adıyaman Province), western Iran and the Iraqi Kurdistan, but samples were available from only three localities (Figure S28).

Helix kazouiniana (Pallary, 1939) was recovered monophyletic, but with two divergent lineages (Figure S29). Even more divergent lineages were found within Helix engaddensis Bourguignat, 1852 (Figure S30), although samples from only two localities just 36 km apart were analysed. Helix engaddensis is a common species in much of Israel and Palestine and western Jordan. Heller [110] distinguished an undescribed form from high elevations at Mount Hermon as a probable separate species (but see [10]) and the darkly coloured forms from east and south of Lake Kinneret as a subspecies Helix engaddensis prasinata Roth, 1855. A similar dark form has been found at the ruins of the crusaders’ castle Krak des Chevaliers in the Homs Governorate in Syria. All these potentially distinct forms are yet to be sampled and assessed phylogenetically.

The remaining two Pelasga species recognized by Neubert [10] are characterized by marked spiral sculpture on the shell: Helix nucula Mousson, 1854 distributed along the south-western coast of Anatolia, some Aegean islands and Cyprus and Helix pachya Bourguignat, 1860 from the Levant. While the former is often small, the latter reaches 5 cm in shell diameter. Surprisingly, the mitochondrial phylogeny showed that these two cannot be separated [23]. There is a clade specific to Cyprus, but samples from the west of Anatolia and Aegean islands on the one hand and from the Levant on the other are intermingled in the tree (Figure S31). The very large form found from Syria to Mt. Hermon in northern Israel (represented in our dataset, however, by just one sampling site in Lebanon) does not seem to be phylogenetically distinct from smaller forms, which are found in the northern Levant, south-western Anatolia and south-eastern Aegean. We therefore conclude that H. pachya is a junior synonym of H. nucula and we suspect that H. nucula is naturally distributed in the Levant, while the distribution in south-western Anatolia and the Aegean may be a result of anthropogenic translocations. Such distribution parallels the cases of Levantina [42] and Helix cincta with its relatives [41]. However, this remains a speculative hypothesis until further detailed phylogeographic analysis is performed in the Levant.

European Clade

The so-called European clade [23] is a well-supported group of Helix species with diversity centre in the western Balkans, comprising eight currently recognized species including the type species of the genus. Relationships between species within the clade remain largely unresolved, but their intraspecific diversity is usually very well sampled. Further detailed sampling may reveal details of contact zones or origin of specific populations. Of the currently recognized species, Helix lutescens Rossmässler, 1837, Helix pomatia Linnaeus, 1758, Helix secernenda Rossmässler, 1847 and Helix thessalica O. Boettger, 1886 are well supported by the data as monophyletic groups and are well defined conchologically. The remaining species, all found in the western Balkans, are more complicated and their relationships cannot be fully resolved with mtDNA data only.

Only shallow divergences were detected within Helix lutescens, which is the only species of the clade distributed exclusively outside the area of the Balkan glacial refugia (Figure S32). It is the only Carpathian biogeographic element within Helicini and, like in H. pomatia or H. thessalica, it apparently performed better in the warmer periods of the last glacial cycles [111,112]. Glacial refugia of H. lutescens may be expected somewhere in Romania, but the existing data do not provide strong hints of their location.

Helix pomatia has been detailed elsewhere [38] and we only provide denser sampling in some areas, especially northern Italy and Croatia (Figure S33). In H. secernenda, our data reveal a centre of diversity in northern Albania and southern Montenegro, in particular around Prokletije Mts. and Lake Skadar, and a late colonization of Dalmatia (Figure S34). For H. thessalica, we provide additional data (mainly from Ukraine) compared to Korábek et al. [40] (Figure S35). Three minor issues remain for this species: the eastward extent of its range is not yet clear, there may be further populations between the southernmost occurrences in Pelion, Greece, and those sampled in Macedonia (see [113]) and the diversity of lineages in a presumed source area for postglacial expansion in the south-western Carpathians [40] needs to be better explored.

A complicated issue is the classification of populations and taxa currently included in Helix dormitoris Kobelt, 1898. So far, this name has been used for several similar forms living predominantly in higher altitudes of the western Balkans [10,114,115,116] and our observations suggest that these may comprise up to four different evolutionary lineages. The typical H. dormitoris is a species from eastern Bosnia and Herzegovina, Montenegro and south-western Serbia (type locality in Durmitor, Montenegro), characterized by a mitochondrial lineage without any close relatives and with a shallowly differentiated crown (Figure S36). Helix dormitoris arnautorum Knipper, 1939 from high altitudes of Šar Planina and Korab yielded a mitochondrial lineage basal to Helix straminea and Helix vladika (Figure S37). Samples from Hajla (Helix dormitoris hajlensis Knipper, 1939) yielded a lineage sister to a small form of H. secernenda from high altitudes of Prokletije (south of Hajla) and is by us provisionally included in H. secernenda (yellow in Figure S34). The last form included in the past within H. dormitoris comes from an isolated limestone massif Mali i Tomorrit in central Albania. The shells there closely resemble typical H. dormitoris or sometimes H. secernenda, but the single tissue sample from the massif yielded a H. schlaeflii haplotype typical for central Albania. The summit area of the massif, where the shells resemble most the other high-altitude “dormitoris” forms, needs to be sampled to resolve the identity of those populations.

Helix straminea Briganti, 1825 and Helix vladika Kobelt, 1898 are very closely related (Figure S37) and there are conchologically intermediate populations. Helix straminea probably colonized the Apennine Peninsula from the territory of current Albania [36]; the present data reveal Albanian populations with an even closer relationship to the Apennine populations than known before. The area along the borders between Albania, Northern Macedonia and Kosovo should be explored in detail to shed light on the degree of isolation between H. straminea, H. vladika and H. dormitoris arnautorum.

Three mtDNA lineages above the “species” threshold exist in Helix schlaeflii Mousson, 1859. Their relationships and the potential monophyly of H. schlaeflii in the mtDNA are unresolved. There do not seem to be readily identifiable conchological differences between individuals of these three clades and their distributions overlap (Figure 5). The northern clade is found in central Albania (Figure S38), roughly north and east of Pogradec, Elbasan and Tiranë. A geographically central clade is broadly distributed in southern Albania, but it has been recorded also from Greece (Figure S39). The third, southern clade was found only in north-western Greece (Epirus, Western Macedonia; Figure S40).

Figure 5. Distribution of “species-level” mitochondrial clades in the European clade of Helix. The distribution of the clades largely corresponds to that of accepted morphospecies except for some samples at range contacts where we assume effects of interspecific hybridization and introgression. The relationships between the three clades found in Helix schlaeflii are unresolved.

Figure 5. Distribution of “species-level” mitochondrial clades in the European clade of Helix. The distribution of the clades largely corresponds to that of accepted morphospecies except for some samples at range contacts where we assume effects of interspecific hybridization and introgression. The relationships between the three clades found in Helix schlaeflii are unresolved.

Our new data on the European clade (Table S1) revealed additional cases of discordance between identifications based on shells and the retrieved mitochondrial lineages to those reported earlier [23,37,38]. Haplotypes belonging to a mtDNA lineage typical for H. pomatia have been found in individuals identified as H. thessalica and H. vladika (Figure S33), lineage of H. vladika in H. dormitoris (Figure S37) and vice versa (Figure S36), lineage of H. thessalica in H. pomatia (Figure S35), lineage of H. dormitoris in H. pomatia (Figure S36), lineage of H. schlaeflii in H. straminea (Figure S38) and lineage of H. secernenda in H. schlaeflii (Figure S34). Sometimes, morphologically intermediate individuals or populations can be found: we have observed these between H. pomatia and H. thessalica, H. pomatia and H. dormitoris, H. vladika and H. dormitoris, H. pomatia and H. vladika. Overall, it seems that interspecific hybridization occurs in very localized contact zones (as we observed in the case of H. pomatia and H. thessalica) because contact of populations of parental species is limited. However, there are areas with populations that are probably of admixed origin as well as possible cases of mitochondrial capture. The discrepancies between shell-based identification and mitochondrial lineages may be a result of contacts between species during range contractions and expansions over the Quaternary glacial cycles [117,118]. This is especially likely in the case of H. dormitoris dormitoris (Figure S36), whose lineage is more broadly distributed than the taxon itself. They were found in H. pomatia and H. vladika at the north-western and southern limits of the range of the H. dormitoris mtDNA lineage, respectively. Helix dormitoris is probably adapted to higher altitudes than these two species and might have been replaced by them after the LGM, leaving behind the introgressed mitogenomes (and probably additional genomic heritage as well, as some of the H. pomatia and H. vladika populations in question have atypical conchological characters). Such events likely occurred also earlier in the evolution of this group, so even if fully resolved, the mitochondrial phylogeny may not fully capture the true relationships between the species.

The European clade of Helix is, besides Codringtonia, the only group where the currently accepted species are distributed in a pattern with a substantial element of parapatric replacement; in other groups the ranges do not adjoin, they overlap more, or the parapatry is limited to two or a few lineages. This is also reflected in the distribution of the “species-level” mitochondrial clades (Figure 5), although their distribution does not fully correspond to that of morphospecies, in particular in the case of the H. dormitoris dormitoris clade. The main exception to the pattern is H. lutescens, which is sympatric with H. pomatia and H. thessalica and is often found in syntopy with these. The frequency of syntopy is not fully apparent from the data because of the difficulty to collect live H. lutescens at some of the visited sites. The ranges of Helix thessalica and H. pomatia overlap mainly in the postglacially colonized areas [38,40] and the two apparently hybridize upon contact [37]. Helix thessalica also lives in sympatry and syntopy with H. vladika, but we saw no phenotypically intermediate individuals. The ranges of Helix straminea and H. schlaeflii overlap in Albania, but no syntopic occurrence is known to us.

The limited overlap of most species’ distributions can be probably partly explained by differing habitat and climate preferences. In particular, Helix pomatia and H. secernenda probably do not come into contact in Croatia as they markedly differ in their climatic niche, the latter being adapted to exposed summer-dry rocky Mediterranean habitats (own observations). Helix schlaeflii has similar preferences. Helix thessalica seems to favour warmer areas than H. pomatia, but both prefer relatively humid sites. There are also forms that seem to be adapted to high altitudes, like H. dormitoris.

Mediterranean Clade

The Mediterranean clade of Helix can be divided into two groups, whose close relationship has only been revealed by molecular phylogenetic analyses: Helix ceratina Shuttleworth, 1843 from Corse and the Apennine complex of lineages related to Helix ligata O. F. Müller, 1774 on the one hand and species related to Helix cincta O. F. Müller, 1774 and Helix melanostoma Draparnaud, 1801 on the other hand [41]. The former group is in a need of a formal taxonomic revision, as Fiorentino et al. [20] only suggested available names for the mitochondrial clades they recovered within the H. ligata complex and not all clades were assigned a name. The authors could not find phenotypic traits that would distinguish members of the different clades. Phylogeny and biogeography of the latter group has been revised recently [41]. The relationships between species are even less resolved than in the European clade, but in this case also the distribution of mitochondrial lineages is poorly known in some species.

Helix ceratina is an extremely threatened species known only from a single small site near Ajaccio, Corse (Figure S41), which is a remnant of a once broader distribution documented by findings of subfossil shells elsewhere on the island [119].

The earliest split within the Apennine Helix ligata complex divides its diversity into a southern and a northern clade (Figure 3). Both are highly diverse. The results suggest marked geographic structuring within the southern clade, but sympatric occurrence of haplotypes from different lineages suggests that this is at least in part due to differences in frequency of individual lineages within populations rather than strict allopatry (Figure S42). Furthermore, the number of sampled populations is too low to reveal the real geographic structure. The distribution of the southern clade overlaps with that of the northern one. Within the latter, lineages overlap in distribution (Figure S43) and we found no apparent geographic structure among them.

Regardless of the systematic status of all mitochondrial lineages, the Helix ligata complex likely contains more than one species. Within the northern clade, there is a peculiar morphotype with white rounded shells with narrow bands, which is, at least in part, associated with high altitudes near or above the treeline (Gran Sasso, Monti Reatini, Majella). It used to be identified as Helix delpretiana Paulucci, 1878 [120,121] and our samples of this form fall within a clade labelled as Helix pomatella Kobelt, 1876 by Fiorentino et al. [20]. In Abruzzo, from where we analysed samples, its range overlaps with the clade considered H. ligata s. str. by Fiorentino et al. [20], but they apparently live in different habitats. However, there appear to be populations of intermediate appearance and the corresponding mitochondrial lineage is more broadly distributed, including lower altitudes. We do not know how closely the distribution of this clade is mirrored by the distribution of the conchological varieties. In addition, another isolated mitochondrial lineage has been recorded from high altitude in Monti Marsicani [20]. We hypothesize that these populations from high altitudes represent relics of lineages more broadly distributed during the Quaternary glacial periods. The spread of H. ligata s. str. then might have led to admixed populations in lower altitudes. Of course, genomic-scale data would be needed to test this hypothesis. Type localities of all nominal taxa in the Apennine Helix need to be sampled to associate phylogenetic clades with available names.

The distribution of several species of the Mediterranean clade has been affected by introductions to new areas within the Mediterranean basin [41]. All but one sequence of H. cincta currently available come from populations outside its native range (Figure S44), which lies probably largely in Syria. In the mitochondrial phylogeny, this species is one of the terminal branches of a “species” level clade from the northern Levant that includes three more recognized species [10]. The two most closely related to H. cincta are not sufficiently known: the typical form of Helix anctostoma von Martens, 1874 from Belen pass in the Hatay province [10] has not been sampled and the known range of Helix valentini Kobelt, 1891 extends to Syria [10], but samples were available only from Turkey. The distribution, monophyly and morphological distinctiveness of these three taxa across their distribution range still need to be established. The combined range of these three overlaps with that of Helix fathallae Nägele, 1901, which is paraphyletic to them in the mitochondrial tree. Its most basal lineage comes from an atypical population and fell above our “species” threshold (Figure S45).

Helix borealis Mousson, 1859 is well sampled. The monophyly of the species is uncertain: it consists of three divergent allopatric clades. The nominotypical form lives in western Greece in three parapatric lineages (Figure S46). The other two clades seem to be relictual and their distribution is limited. One is found on northern Evvia and in Northern Sporades (Figure S47), the other in two isolated areas on Crete and in south-western Anatolia (Figure S48).

The two African species, Helix melanostoma and Helix pronuba Westerlund & Blanc, 1979, are also poorly sampled. For H. melanostoma (Figure S49), there is no sequence available from Algeria, which comprises approximately two thirds of its range in the west. The data for H. pronuba (Figure S50) are fragmentary, as no fresh or properly preserved samples were analysed. The core of its Holocene broad range may be Cyrenaica, with subfossil shells dated to 17,000–14,000 cal. BP [122].

Anatolian Clade

The Anatolian clade comprises six currently accepted species [23]. All live in Anatolia, although the ranges of two extend also into south-eastern Europe and one lives also in the Caucasus.

Helix asemnis Bourguignat, 1860 is, compared to other Helix species, currently broadly delimited in respect to conchological variation and the divergences between its mtDNA lineages [10,23]. It consists of two clades, eastern (Figure S51) and western (Figure S52), whose ranges adjoin along the Ecemiş fault zone in southern Turkey (roughly along the line Pozantı–Mersin; a corresponding divide probably exists also in Levantina cilicica (Kobelt, 1895), Figure S11). Both clades are diverse and the individual lineages within both seem distributed in an allopatric manner. There is also a substantial corresponding geographic conchological variation, suggesting that H. asemnis may comprise several narrowly distributed species. As in the H. ligata complex and Codringtonia, finer sampling would be needed to reveal the degree of isolation or distribution overlap between the uncovered lineages. We suggest that the north-eastern and altitudinal range limits of H. asemnis should also be better explored.

A species closely related to H. asemnis is Helix escherichi O. Boettger, 1898 from north-western Anatolia (Figure S53). Its current range is only poorly known and there are no recent samples available from anywhere near its type locality Akşehir in the Konya Province [10]. As regards Helix pathetica Mousson, 1854, the species has a very broad distribution range in inner Anatolia, which is probably young, perhaps of only Holocene age, because of the very small differences between haplotypes from across central Anatolia (Figure S54). We hypothesize that more mitochondrial diversity within the species might be found in the as yet unsampled north-west of its range, because north-western Anatolia appears to be a diversification centre of the Anatolian clade. Anatolian samples are missing altogether for Helix pomacella Mousson, 1854 (Figure S55), which lives around the Sea of Marmara and up to Burgas [10] and is so far represented by only one analysed individual from Bulgaria.

The divergence between the mitochondrial clades of Helix nicaeensis Férussac, 1821 and Helix lucorum Linnaeus, 1758 is in a range observed within some species (Figure S56). However, the two species are morphologically so different that H. nicaeensis has even been placed into its own genus Tacheopsis prior to molecular phylogenetic analyses [123]. Helix lucorum has been recently extending its range through introductions to anthropogenic habitats ([39] and references therein) and data posted online at the iNaturalist website (https://www.inaturalist.org (accessed on 25 December 2021)) suggest a far greater extent of the introductions than covered by peer-reviewed literature, for example in Central Asia. The present data on the distribution of its intraspecific lineages are nevertheless insufficient for identifying the geographic origins of the lineages involved in the expansion. Published figures and photos posted on iNaturalist show that most of the newly emerging non-native populations are of a morphotype distributed from Europe along the southern Black Sea coast up to the western Caucasus and associated with a specific mitochondrial lineage [39]. The natural distribution limits of this lineage Europe and western Caucasus are disputed [39,124], leaving northern Anatolia as its possible cradle.

Species Unassigned to Clades

Four Helix species cannot be assigned to any of the above four major mitochondrial clades. Helix buchii Dubois de Montpéreux, 1840 lives in north-eastern Turkey, in Georgia and in part of Armenia [10]. Like in H. nucula, there are remarkable differences in shell shape and size between populations (3–6 cm in diameter [10]): individuals from some of those located westerly are very small (e.g., from Espiye) while other populations more to the east make the species the largest helicid. Two divergent, unrelated lineages have been recovered from H. buchii. One lineage, belonging to the Anatolian clade but distinct from other species, was recovered from a single individual collected near the Sümela monastery in the Trabzon Province of Turkey (Figure S57). All other H. buchii individuals analysed so far yielded an unrelated lineage outside the Anatolian clade, but given that the deviating individual shared the shell characteristics typical for H. buchii in the same region, we consider unlikely that the sample represents a distinct species. This case may represent a “ghost” mitochondrial lineage from a past introgression, the source of which is either extinct or yet to be uncovered. The data from the remaining samples point to a decrease in diversity from west to east (Figure S58) with an overlapping distribution of mitochondrial lineages, but remain too scarce for a reliable description of the phylogeographic pattern.

Helix pelagonesica (Rolle, 1898) has a small range extending from Thessaly and Macedonia, Greece, to south-east North Macedonia. The available data indicate substantial mtDNA diversity but are insufficient for description of potential geographic structure within its small range (Figure S59). However, an isolated locality lying outside the known range is reported here from Morfovouni near Karditsa in Central Greece, which yielded a lineage basal to the other H. pelagonesica samples.

The current range of Helix philibinensis Rossmässler, 1939, stretching from Lake Prespa in the west to Asenovgrad near Plovdiv in Bulgaria is apparently young, as samples from all the range extremes yielded identical or very similar haplotypes (Figure S60). Interestingly, its distribution appears to be in large part patchy, which cannot be explained by patchiness of suitable habitats. We have found H. philibinensis on different bedrock from limestone to granite and although it prefers warm open habitats, we found it also in shaded places under tree cover.

Helix albescens Rossmässler, 1839 has a large range stretching from Azerbaijan to Ukraine (Mykolaiv, Odessa), but sequence data from autochthonous populations are available only from the Caucasus and Crimea (Figure S61). The single presented sample from Bulgaria (an individual from a site without GPS coordinates near Ivaylovgrad [10,23]) is problematic. Helix albescens is not known anywhere else in Bulgaria (I. Dedov, pers. comm.) or from Romania [125]. Neubert [10] lists additional Bulgarian localities, but we did not find H. albescens at any of these and, upon inspection of the original material, we conclude that these records are based on misidentified material. We did not find the species in the immediate vicinity of Ivaylovgrad (instead, H. figulina was present), but the sampling site probably lies farther from the town by the Ivyalovgrad Reservoir (Table S1).

3.5. Distribution of Intraspecific Diversity

Even though insufficient data are available for many species at the moment, several observations can already be made regarding the intraspecific diversity and its distribution that emerge from the comparison across the tribe (Figure 6).

Clades of similar age vary profoundly in their range size; largest range dimension between small and large ranges may differ by two orders of magnitude. The largest ranges of intraspecific lineages are mostly located in central and eastern Europe and result from postglacial range expansions [38,40]. A specific case is the only broadly distributed intraspecific lineage within H. lucorum, the spread of which was apparently greatly facilitated by anthropogenic dispersal even before the documented establishment of many newly founded populations since the late 19th century [39]. This species has been collected and transported for food, but its spread is clearly facilitated also by its broad ecological valence and tolerance for anthropogenic habitats. Similar ecological characteristics helped the expansion of H. pomatia. Tolerance to soils with relatively low calcium content might be an advantageous trait for both species. A factor contributing to expansion of H. thessalica and C. vindobonensis could be also the relative topographic homogeneity and zonal distribution of biomes in eastern Europe in comparison with the geographically complex Balkans.

In contrast to the large ranges of the postglacial colonizers, there are whole species with extremely restricted distributions (A. praecellens, N. neocrassa, L. christophi). Amanica praecellens and N. neocrassa are particularly isolated old lineages (Figure 3), whose ranges may be relictual. There are also somewhat more broadly distributed taxa with strong intraspecific geographic structure of narrowly distributed mitochondrial lineages (H. asemnis, Codringtonia). Both situations require long-term environmental stability allowing for differentiation and persistence of the lineages in question and are thus informative about the history of the respective regions. There are several areas with a pattern of allopatric, narrowly distributed lineages (or where this pattern may be suspected based on the current data). One such region is in the Taurus Mts. in Cilicia, southern Turkey, where similar structure is observed in sympatric H. asemnis and L. cilicica. Additional such areas are found in Europe. These include the Peloponnese (with Codringtonia), Albania (H. schlaeflii, H. straminea, H. secernenda), southern Apennine Peninsula (H. ligata complex) and likely the western edge of the Greater Caucasus (C. atrolabiata).

It is not uncommon that two or more divergent intraspecific lineages are distributed sympatrically and even co-occur within local populations. This may happen in postglacially colonized areas thanks to colonization from multiple sources (H. pomatia; [38]), but the same pattern may be found also in diversifications centres where the species survived through the glacials. The pattern in H. pomatia combines a geographic structure observed in Bosnia and co-occurrence of several lineages in potential refugia in western Romania [38]. Additional possible examples of high diversity in a stable range core are not yet sufficiently sampled (H. pronuba, C. leucoranea). Geographic structuring means lower local lineage diversity than when the lineages are sympatric, but both patterns may indicate regional environmental stability; the situations differ in the population connectivity across the region. The cases of the finest geographic structure like in H. asemnis or Codringtonia would require extremely dense sampling to evaluate whether boundaries between the ranges of the intraspecific clades are sharp or rather overlapping. In addition, the sample sizes per sampling site would have to be larger to capture the rare lineages within each sampled population, because what appears to be a strict allopatry may in fact be just differences in lineage frequencies. These characteristics are relevant for answering questions regarding the role of dispersal barriers in the origin and persistence of phylogeographic structures.

An important point to consider is the range of spatial dynamics suggested by the available data. On the one hand, there are species comprising deeply divergent lineages which might have persisted in the same region for millions of years (possible examples: M. dickhauti, H. kazouiniana, H. engaddensis); on the other hand, there are some remarkable postglacial expansions (H. pomatia, H. thessalica, C. vindobonensis) and some relatively old species that have recent crown ages (C. vindobonensis, H. philibinensis, H. lutescens, possibly H. pathetica). While the stem ages of these species are around 10 My or more, their crown groups may be more than ten times younger than their stems (depending on the magnitude of time-dependency of the clock rate [126,127]). If not caused by selection, this could be a result of past range contractions, as documented for the glacial cycles in central and Eastern Europe. For example, fossils attributable to the stem lineage of C. vindobonensis occur in deposits of the late Miocene age (most likely Tortonian [128]) near Kavarna in eastern Bulgaria (own observations) and its closest relatives live in the Caucasus and Alborz [24], yet the crown group likely originated in the western Balkans during the second half of the Pleistocene and large parts of the current distribution range of C. vindobonensis date to less than 12,000 years [40]. That indicates that most of the biogeographic history of these species captured in mitochondrial genomes has been lost due to extinctions and this hidden past might have involved substantial and repeated changes in range extent and position.

3.6. Outlook

The phylogenetic and phylogeographic studies based on sequences of mitochondrial markers revolutionized the understanding of (not only) land snail diversity and its roots and had a tremendous impact on taxonomy. That holds also for the group detailed here and it is very likely that for some time the mitochondrial data will remain indispensable.

3.6.1. Mitochondrial Data in Helicini Taxonomy

Despite some persisting problems, the potential of the mitochondrial-only data for taxonomy and phylogenetics is now nearly exhausted in Helicini. They enabled great progress in the systematics of Helicini, allowing to sort the many described forms (e.g., [10]) into natural groups. Genera are now reliably delimited, groupings of closely related species revealed, several species redefined and taxa warranting further systematic attention identified. However, the limitations of the mitochondrial data are obvious (see [129] for an illustrative example) and we expect that with detailed studies using genomic-scale data and considering ecology of the snails, the species-level taxonomy would see additional changes.

There are subtle indications that the biological species are in some cases smaller units than recognized by the current taxonomy. While the intraspecific mitochondrial diversity may be substantial and old (e.g., in Helix pomatia, H. thessalica), there may also be species characterized by young mitochondrial lineages and some of them may have indeed formed recently. In central Italy, there appear to be specialized high-latitude forms which may deserve recognition as distinct species. In the Balkans, the high-latitude populations from Prokletije and Hajla, placed here in H. secernenda, probably also have a different set of adaptations than the typical H. secernenda from the warm Dalmatian coastland. A particular case presents H. asemnis, where the locally distributed intraspecific mitochondrial lineages seem to be associated with specific morphotypes, indicating differentiation far beyond the mitochondrial genomes. Another similar case may be H. cincta with its relatives, but the distribution of the conchological forms known as H. cincta, H. valentini and H. anctostoma and their association with a particular lineage needs to be clarified with good sampling from the Hatay region and Syria. Finally, crossing experiments showed reduced fitness of hybrids between populations classified as subspecies of L. spiriplana [130], also suggesting an advanced stage of speciation.

3.6.2. Mitochondrial Data in Phylogeography

Further sequencing of selected mitochondrial genes remains the most feasible approach to learn more about the variation in phylogeographic histories and about the distribution of diversification centres and refugia of land snails in the Western Palaearctic. Taxonomic and geographic coverage offered by the mitochondrial data in land snails is not in sight with multilocus data, for reasons that include not only the costs of alternative methods (decreasing but still substantial) and their complexity, but especially the availability of suitable samples. The mitochondrial data accumulated over many years of research thus represent a unique resource, which would be reasonable to improve further by expanding the geographic and taxonomic coverage in understudied regions to fully exploit its potential. Furthermore, taxonomy, to serve its purpose, must maintain its continuity dating back to the second half of the 18th century and new taxonomic research should relate to earlier hypotheses. The data compiled here include material used for earlier morphology-based taxonomic revisions of the genera Codringtonia, Levantina, Isaurica and Helix [10,95,96,101], allowing for connecting taxonomic hypotheses with phylogenetic lineages, though only mitochondrial ones. This makes it meaningful to include an analysis of mitochondrial sequences also in parallel to future genomic-scale studies in order to provide the link from the new types of data to the earlier work. We recommend that both most commonly used fragments of mitochondrial genes for 16S and cox1, defined by the primer pairs 16Scs1+16Scs2 and LCO1490+HC02198 (Table 1), respectively, are targeted in such cases, with priority given to 16S when only one is used. However, adding additional loci or sequencing of whole mitogenomes would be helpful for some research questions. Despite the wealth of data, well resolved mitochondrial phylogeny is lacking, preventing for example comparison of mitochondrial and nuclear phylogenies and divergence dates. For phylogeography, including faster evolving mitochondrial genes could be helpful as the variability in the relatively slowly evolving 16S and cox1 limits their use for locating glacial refugia.

The sampling density needed to uncover the regional diversity and its structure widely differs between regions, but for an unbiased view of the distribution of intraspecific diversity, it is vital that species and areas are sampled indiscriminately at least in the beginning, without preference for species with problematic taxonomies or where readily interpretable phylogeographic structure is expected. The data collated here show a broad range of different phylogeographic structures and histories. They suggest several relatively recent range expansions of different extent and uncover several regional diversification centres or refugia. If the sampling focused only on well-established refugia and diversity centres, the full extent of variability in population histories and the dynamic nature of the distribution of many of the species and intraspecific lineages would be obscured. However, very fine sampling is still necessary to characterize the distribution of lineages in the diversity centres, as in the case of H. asemnis or the H. ligata complex.

Further progress in documenting the diversity of Helicini heavily depends on collecting new samples from the eastern half of the tribe’s distribution. That holds in particular for Isaurica and Levantina, where even some of the currently accepted species are not sampled at all and for all Helicini members in the Levant, eastern Turkey and adjacent regions of Iraq and Iran. It is difficult, or almost impossible, to move forward without a broad participation of local zoologists. Geographically complete sampling is easiest for locals, who may also leverage on their knowledge of the regional biogeography to identify populations worth sampling (e.g., by considering known regional dispersal barriers and diversity hotspots) as well as phenology, distribution of suitable habitats (Figure S62) and other factors when planning the sampling. We are open for cooperation and willing to help anyone interested in the diversity of the Helicini.