Diversification over deep and shallow temporal scales in the Holarctic genus Perpolita (Gastropoda: Gastrodontidae)

Abstract

The Holarctic land snail genus Perpolita was used to explore the influence of past and current biogeography on diversification. The number of empirically-supported species was determined using a consensus between mtDNA sequence, nDNA sequence, conchology, and geographic and ecological range with five valid temperate-boreal species (Perpolita binneyana, Perpolita electrina, Perpolita hammonis, Perpolita petronella, and Perpolita radiatella) being recognized. Only P. petronella was unchanged in both nomenclature and diagnostic characteristics with the remainder requiring alterations. Perhaps the most important of these was elevation of P. radiatella to a valid species, with its populations having been previously lumped either under European P. hammonis or North American P. electrina. Divergence times of 18.7–10.0 Mya were suggested through genome-wide SNPs in combination with the fossil record, indicating a pre-Pleistocene origin for all Perpolita species. Using genetically-confirmed diagnostic shell characters, we accumulated > 2000 valid occurrences and used these to estimate appropriate modern and Last Glacial Maximum climate extents for all species. These models suggest that modern intra-specific gene pool diversity may generally reflect Pleistocene palaeoclimatology.

biogeography, ddRAD-seq, integrative revision, Last Glacial Maximum, morphology, Nesovitrea, taxonomic error

Introduction

Current biodiversity is not only related to modern environmental gradients, ecological interactions, and biogeography but also to the historical template (e.g. Ricklefs and Schluter 1993, Nekola and White 1999). While cyclic Pleistocene climate change has been implicated in diversification (e.g. Weir and Schluter 2004, Hope et al. 2012, Weir et al. 2016), many species’ origins pre-date the Pleistocene (e.g. Toews and Irwin 2008, Kleckova et al. 2015, Schär et al. 2018). The diversification process must thus consider a mix of both ancient and recent drivers.

The land snail genus Perpolita H.B. Baker, 1928 appears potentially useful in considering such issues. It likely originated in the Oligocene (Schlickum and Strauch 1975), and is currently believed to represent six nominal species and one subspecies across the Holarctic (MolluscaBase 2024). While some authors have placed this genus within Hawaiian Nesovitrea C.M. Cooke, 1921 (e.g. Hubricht 1985, Sysoev and Schileyko 2009, Welter-Schultes 2012), this act has not been based on any empirical tests. Hausdorf (1998), Bouchet et al. (2017), and MolluscaBase (2024) keep these genera separate, and we follow that approach here. Two species [Perpolita petronella (L. Pfeiffer, 1853), Perpolita hammonis (Strøm, 1765)] are thought limited to Eurasia (Sysoev and Schileyko 2009, Welter-Schultes 2012), two [Perpolita binneyana (E.S. Morse, 1864); Perpolita electrina (Gould, 1841)] to temperate/boreal North America (Pilsbry 1948, Hubricht 1985), with the remaining two [Perpolita dalliana (Pilsbry & Simpson, 1889), Perpolita suzannae (Pratt, 1978)] being limited to warm temperate/subtropical areas adjacent to the Gulf of Mexico in southeastern North America (Pilsbry 1948, Hubricht 1985). In addition, P. binneyana has been segregated into two subspecies with Perpolita binneyana occidentalis H.B. Baker, 1930 being used to demarcate populations west of 100°W (Pilsbry 1948).

This genus demonstrates parallel shell traits in species restricted to different continents, with P. electrina/P. hammonis both possessing darker shells and P. binneyana/P. petronella possessing lighter/uncoloured shells. Additionally, while the ecological niches for all four boreal taxa overlap to some extent, P. binneyana/P. hammonis both tend to occupy drier sites than P. electrina/P. petronella (Nekola 2004, Welter-Schultes 2012, Horsák et al. 2013). Given that empirical testing of these taxonomic hypotheses has never been attempted, it is unclear whether the patterns are real and require an evolutionary explanation, or whether they simply reflect unjustified taxonomist exuberance within phenotypically or ecologically plastic taxa (e.g. Nekola and Horsák 2022).

We investigate these issues through a global analysis of all currently recognized Perpolita taxa. We first determine the number of empirically supported taxa using mtDNA sequence, nDNA sequence, and shell morphology as our distinct data streams. We then use nDNA SNP data in combination with dated fossils to document species-scale diversification. We use these genetically-verified concepts to determine robust shell identification features, and then use these characters to generate valid range and occurrence information for each species. These data were used to estimate the geographic extent of modern and Last Glacial Maximum (LGM) climate niches of each species through ecological niche modelling (ENM). The impact of past and present ranges on the modern intraspecific gene pool diversity was then considered.

Material and methods

Taxon sampling and DNA methods

We analysed 111 individuals representing the entire ecological and biogeographic range of all currently recognized Perpolita taxa. These were sourced from 22 countries and generally collected from 1995–2022 (Supporting Information, Table S1; Fig. 1). Specimen selection was undertaken through an iterative process: we initially assembled a set of 30 individuals (e.g. 3–6 per taxon) across the known geographic and ecological range of all taxa using currently accepted taxonomic concepts and diagnostic characters (Pilsbry 1948, Sysoev and Schileyko 2009, Welter-Schultes 2012, Hayase et al. 2016). These were subjected to an initial phylogenetic analysis, with concordance/discordance between prior and actual highly-supported concepts being noted. We then sorted material into genetically-validated groups and used these to determine potential robust diagnostic shell characters. The ability of these to accurately demarcate genetic clades was then tested by using them to select additional individuals from other populations. These then had their identifications verified using both mtDNA and nDNA sequence barcodes. This procedure was repeated until no mismatch occurred between morphological and genetic-based identifications. Documented robust traits were then used to select a final group of individuals to ensure the total specimen set included the full biogeographic and ecological extent of each genetically confirmed entity. Following selection, each specimen was imaged (see below) with total DNAs then being extracted using the E.Z.N.A. Mollusc DNA Kit (Omega BioTek, GA, USA) following manufacturer protocols.

Figure 1.

Distribution of Perpolita occurrences used in this study. The black crosses indicate all examined samples, and the coloured dots represent samples used in the molecular analyses. Map retrieved from Natural Earth (https://www.naturalearthdata.com/), and presented with the azimuthal equidistance projection centred on the North Pole using QGIS 3.16 (QGIS Development Team 2020).

We amplified and sequenced one mitochondrial [cytochrome b (Cytb)] and two nuclear regions [Internal Transcribed Spacer 1 of the rRNA gene cassette (ITS1) and Intron 8 of the Embryonic Lethality and Abnormal Visual System (ELAV8; Nekola et al. 2022b)] per specimen. Primers and amplification protocols are provided in Table 1. Following ExoSAP purification of PCR products, Sanger sequencing was conducted on each by the SEQme s.r.o. laboratory (Czech Republic) using the Gerbera Sequencing Kit v.3.1. The primer ends were removed, and all resulting sequences were deposited in GenBank (accession nos PP565243–PP565334 and PP573309–PP573503).

Table 1.

Information on primers and PCR conditions used in this study. See also Nekola et al. (2022b) for PCR of the ELAV8 region

Primer	Direction	Sequence 5’–3’	PCR conditions	Reference
Cytb
CytB_Nes_fe	Forward	GGA CGT GGA ATT TAT TAT CAA AG	96 °C 10 min, (94 °C 30 s, 45 °C 30 s, 72 °C 60 s) × 40, 72 °C 10 min	Authors
CytB_NesUni_r	Reverse	GCR AAT AAR AAA TAT CAC TCR GG		Authors
ITS1
ITS1-F	Forward	TAA CAA GGT TTC CGT ATG TGA A	96 °C 10 min, (94 °C 30 s, 52 °C 30 s, 72 °C 60 s) × 40, 72 °C 10 min	Armbruster et al. 2000
ITS1-R	Reverse	TCA CAT TAA TTC TCG CAG CTA G		Horsáková et al. 2019
ELAV8
ELAVF_NI_100	Forward	ATG ACA CAG TTG GAC TTG GAG	96 °C 10 min, (94 °C 30 s, 52 °C 30 s, 72 °C 60 s) × 40, 72 °C 10 min (slightly modified by authors)	Nekola et al. 2022b
ELAVF_NI_84	Reverse	GAG AAC AAC TTC TCC AAA TCC T		Nekola et al. 2022b
ELAVF_NO_88	Forward	GCT AAC TTA TAT ATT AGT GGC TTG C	96 °C 10 min, (94 °C 30 s, 50 °C 30 s, 72 °C 60 s) × 40, 72 °C 10 min	Nekola et al. 2022b

Molecular phylogenetics

Phylogenetic analyses were independently performed on each locus to detect potential topological incongruences using MEGA 7 (Kumar et al. 2016). MUSCLE (Edgar 2004) was employed to align each dataset with default parameters. Final alignments were adjusted by eye. To reduce phylogenetic uncertainty, alignment gaps in ELAV8 and ITS1 were eliminated using trimAl 1.2 (Capella-Gutiérrez et al. 2009). Identical haplotypes were collapsed with FaBox 1.61 (Villesen 2009). Maximum likelihood analyses (ML) were conducted using IQ-TREE v.2.2.0 (Minh et al. 2020). PartitionFinder 2 (Lanfear et al. 2017) was employed to select evolutionary models and partitions based on corrected Akaike information criterion (AICc; Table 2). Ultrafast bootstrapping (Hoang et al. 2018) with 1000 replicates was utilized to evaluate the ML tree topologies. Bayesian inference analyses (BI) were performed using MrBayes 3.2.7a (Ronquist and Huelsenbeck 2003), with four simultaneous chains and sampling trees every 1000 generations. All datasets were run for 1 million generations with four chains. The first 101 trees were discarded as burn-in after assessing convergence with Tracer v.1.7 (Rambaut et al. 2018) and ensuring that the effective sample size (ESS) of all parameters in each dataset was greater than 200. The remaining trees were then used to estimate each maximum clade credibility tree and to evaluate each tree topology with Bayesian posterior probability (BPP) derived from the posterior samples using MrBayes 3.2.7a.

Table 2.

Evolutionary model and partitioning for each phylogenetic analysis based on PartitionFinder 2 estimation (Lanfear et al. 2017)

Alignment	Evolutionary model and partitioning for IQ-TREE	Evolutionary model and partitioning for MrBayes
Cytb	Cytb_pos1: TVM+G Cytb_pos2: GTR+G Cytb_pos3: F81	Cytb_pos1: GTR+G+X Cytb_pos2: GTR+G+X Cytb_pos3: HKY+X
ITS1	TVM+I+G	GTR+I+G
ELAV8	TVM+I	GTR+I
ITS1+ELAV8	ITS1: TVM+I+G ELAV8: TVM+I+G	ITS1: GTR+I+G ELAV8: GTR+I+G
ddRAD loci	K3Pu+F+I+G4	-

Because ITS1 and ELAV8 showed no topological incongruencies within highly supported branches (Supporting Information, Figs S1–S3), both were concatenated into a single nDNA construct on which the following phylogenetic analyses were based. However, because of topological inconsistencies, mtDNA and nDNA data streams were analysed separately with consensus patterns being identified.

Genome-wide Single Nucleotide Polymorphisms methods

We also performed genome-wide Single Nucleotide Polymorphism (SNPs) analysis to estimate phylogeny across a broader part of the nuclear genome. This was limited to 16 individuals spread across the five empirically validated Holarctic species-level groups identified in the above analyses with subtropical P. dalliana being used as the outgroup (Supporting Information, Table S1). Genome-wide SNPs were obtained through double digest restriction-site associated DNA (ddRAD) sequencing (Peterson et al. 2012). This method has recently been used to study a variety of non-model organisms including molluscs (e.g. Haponski et al. 2017, Sano et al. 2022, Hirano et al. 2023). To construct the sequencing library, a double digest was performed using two restriction enzymes (EcoRI and MspI) with approximately 50 ng of extracted genomic DNA. P1 and P2 adapters from Peterson et al. (2012) were ligated to the digested DNA fragments. After the ligation, the DNA fragments were multiplexed and then purified using NucleoSpin Gel and PCR Clean-up kit (Macherey-Nagel, Germany). The library was then amplified for 10 PCR cycles. PCR conditions followed the protocol of Peterson et al. (2012), with final primer concentration being modified to 0.5 µM. Finally, the amplified sample was purified using the same clean-up kit. DNA fragments of 200–400 bp were collected using Pippin Prep (Sage Science, MA, USA). The constructed DNA library was sent to Macrogen Japan (Tokyo) and sequenced using an Illumina Hi-seq X paired-end 150 bp setting (Illumina, CA, USA). The sequenced raw reads were then processed for de novo assembly using ipyrad v.0.7.29 (Eaton and Overcast 2020). We specified a minimum coverage of nine (> 50%) samples per locus, with strict filter adapters (filter_adapters = 2) to completely remove the Illumina adaptor sequences. Other parameters were set to default which are adjusted by the developers in most cases. We were able to analyse 78 717 SNPs distributed across 1 171 157 bp in 4134 loci. The raw sequences were deposited in the NCBI Sequence Read Archive (PRJNA1094558).

The obtained SNPs were then used to conduct two different phylogenetic analyses: first, ML reconstruction was based on concatenated sequences using IQ-TREE v.2.2.0 (Minh et al. 2020). The nucleotide substitution model was automatically selected using ModelFinder Plus (Kalyaanamoorthy et al. 2017) as implemented in IQ-TREE (Table 2). ML tree topologies were evaluated using ultrafast bootstrapping (Hoang et al. 2018) with 1000 replicates. Second, SVDquartets (Chifman and Kubatko 2014) was implemented in PAUP * v.4.0a168 (Swofford 2002). This method estimates a tree using SNPs data under the coalescent model, unlike analysis from concatenated sequences, which may erroneously lead to high support values (Gadagkar et al. 2005), and furthermore, unlike methods using gene trees for each locus, has high statistical consistency and is considered robust against discrepancies between gene trees (Chifman and Kubatko 2014, Wascher and Kubatko 2021). All possible quartets were sampled, and the quartet Fiduccia and Mattheyses algorithm was employed for tree inference (Reaz et al. 2014). Bootstrap sampling (BS) with 1000 replications was performed to test the topology.

Divergence time estimation

We used Perpolita wenzi Schlickum & Strauch, 1975 from the lowermost Miocene Aquitanian period (23–20.4 Mya) in Erbach, Germany, as the oldest known Perpolita (Schlickum and Strauch 1975). This entity has clear Perpolita characters such as regular transverse striae on the upper surface, giving the shell a radiated appearance on gradually expanding coiling. We excluded the nominal oldest fossil—Perpolita mattiaca (O. Boettger, 1903) from the uppermost Oligocene Chattian period (23–28 Mya)—because of its wide umbilicus (1/4 of the shell width), obtusely carinated shell with lack of transverse striae, and high coil expansion rate (Boettger 1903). These characters define the shell morphology of the genus Aegopinella with the small size and carinated periphery making it similar to modern Aegopinella pura (Alder, 1830). In fact, A. pura was used as the comparative taxa in the original P. mattiaca description (Boettger 1903). The calibration was set as a lower constraint on the node of the most recent common ancestor (tMRCA) of the Holarctic Perpolita. The uncertainty of the upper constraint node (i.e. tMRCA for all Perpolita) was the reason why the lower constraint was used in the analysis. This uncertainty is due to distinct shell morphology (small size and carinated periphery) of P. dalliana/suzannae compared to P. wenzi. Therefore, we assigned P. wenzi to the branch before the node of the tMRCA node of Holarctic Perpolita. Because the phylogenetic relationship between P. dalliana/suzannae and the other Perpolita species remains uncertain (see below), the calibration was set as a lower constraint on the crown node of Holarctic Perpolita. We utilized a lognormal distribution with 23 Mya as the upper limit of the 95% confidence interval (CI; 23.00–16.03 Mya).

We employed two different methods to generate time-correlated trees. First we used RelTime (Tao et al. 2020) which represents a relaxed clock approach using branch lengths. This provides for computationally effective time estimation across all nodes (Tamura et al. 2012, 2018). RelTime was implemented in MEGA X (Stecher et al. 2020) using the above 95% CI and the GTR+I+G substitution model.

Second, we also used BPP 4.4.0 (Rannala and Yang 2003, Flouri et al. 2018). This estimates divergence time under a multi-species coalescent model without introgression using all SNPs loci and assuming GTR base pair substitution. A relaxed correlated molecular clock (Flouri et al. 2022) was employed, with a prior distribution of theta = inverse gamma (3, 0.002) and root tau = inverse gamma (3, 0.08) being adopted. Two hundred thousand Markov chain Monte Carlo samples were collected every two iterations after the initial 10 000 burn-in. Tracer v.1.7 confirmed that the ESS for all parameters was greater than 100. Finally, absolute divergence time with 95% CI for each node was calculated using the R package ‘bppr’ (Angelis and Dos Reis 2015, R Core Team 2022).

Comparative conchology

We recorded maximum shell width and height, umbilicus width, maximum numbers of whorls, shell colour, and strength of shell bottom grooves from 203 individuals sourced from 76 genetically validated populations across all valid species (Supporting Information, Table S1). Shells were imaged using an Olympus szx-7 microscope with an Olympus C-7070 Wide Zoom camera and QuickPHOTO MICRO software at c. 20× magnification. Images were depth composited using Combine ZM (https://combinezm.en.lo4d.com/windows). Microsculpture on lower shell surfaces was imaged using a Keyence vhx-5000 digital microscope with zs-20 and zs-200 objective lenses at 200× magnification.

Climatic niche modelling

Using robust diagnostic features generated through the integrative empirical revision (see above), species identification was verified for all 1781 Perpolita lots in the author’s collections. As P. radiatella is absent from the United Kingdom, we assume that all P. hammonis records provided by Kerney (1999) to be correctly assigned. However, the co-occurrence of P. hammonis and P. radiatella in Sweden prevented the use of P. ‘hammonis’ records from Waldén (2007).

These validated modern observations were then used to parameterize ENMs using protocols described in Nekola et al. (2022a). Climate data at 2.5-minute resolution was retrieved for each occurrence using WorldClim v.2 (Fick and Hijmans 2017). To avoid excessive multicollinearity, climate data was limited to seven variables relevant for land snail distribution (Bio3 = isothermality; Bio4 = temperature seasonality; Bio5 = maximum temperature of warmest month; Bio6 = minimum temperature of coldest month; Bio12 = annual precipitation; Bio15 = precipitation seasonality; Bio18 = precipitation of warmest quarter). Records were pruned to remove redundant observations until the smallest average pairwise distance in reduced PCA space was 0.05 standard deviations. The following data were used for modelling (total recorded/pruned records): P. binneyana = 452/180, P. electrina = 453/173, P. petronella = 280/135, P. hammonis = 348/203, and P. radiatella = 248/147 (detailed coordinates provided in Supporting Information, Table S2).

Climatic niches were generated for each species based on the protocols outlined in Nekola et al. (2022a), using the Maxent algorithm implemented in ‘maxent.jar’ (Hijmans et al. 2021) for the algorithm of ‘ENMevaluate’ function on ENMeval 2.0 (Kass et al. 2021). Maxent is a well-developed method for preforming ENM without absence data (Phillips et al. 2006, 2017) and is widely used in species distribution modelling, with practical instructions (e.g. Warren and Seifert 2011, Merow et al. 2013). It is generally considered to have high accuracy (Kaky et al. 2020, Perkins-Taylor and Frey 2020, Ahmadi et al. 2023) and is particularly robust to sampling bias, which is common in field data (Ahmadi et al. 2023). Random 4-fold partitioning was used for cross-validation. As there is no absolute criterion for model selection we used area under the curve values (AUC; Lobo et al. 2008) and AICc (Warren and Seifert 2011). The best-fitting model for each species (Table 3) was then used to estimate potential range in current and LGM landscapes using dismo 1.3-9 in R (Hijmans et al. 2021). LGM climate data was predicted based on Community Climate System Model 4 (CCSM4; Gent et al. 2011), MIROC-ESM (MR; Watanabe et al. 2011), and MPI-ESM-P (ME; Jungclaus et al. 2011), and were retrieved from WorldClim v.1.4 (Hijmans et al. 2005).

Table 3.

Selected best models for each species in ecological niche modelling. The meaning of each letter in the feature class (FC) is as follows: L = linear, Q = quadratic, P = product, T = threshold. Regularization multipliers (RM) from one to eight were considered. See also Supporting Information (Table S4) for full information on all models

Species	Best AUC model	Average AUC	Best AICc model	AICc
P. binneyana	FC: LQPT/RM: 1	0.907543365	FC: LPT/RM: 2	4800.67003
P. electrina	FC: LQPT/RM: 1	0.933310521	FC: LQT/RM: 2	4432.596038
P. radiatella	FC: LQP/RM: 1	0.881006354	FC: LQPT/RM: 2	4130.013133
P. hammonis	FC: PT/RM: 1	0.96346953	FC: PT/RM: 2	5014.046584
P. petronella	FC: QPT/RM: 1	0.924949036	FC: LQT/RM: 2	3637.574817

Results

Molecular phylogenetics

With full support all phylogenies showed that North American Gulf Coast P. dalliana and P. suzannae are highly distinct from the remaining higher-latitude Perpolita species (Figs 2–3; Supporting Information, Figs S1–S3). However, they possess very small intra-taxon divergence in both their mtDNA and nDNA sequence. Because they appear to represent a distantly-related evolutionary linage, we limit the bulk of the following analyses to the remaining temperate-boreal Perpolita taxa.

Figure 2.

Maximum likelihood phylogenies of the genus Perpolita inferred from mitochondrial cytochrome b, and two nuclear regions [i.e. Internal Transcribed Spacer 1 of the rRNA gene cassette (ITS1) and Intron 8 of the Embryonic Lethality and Abnormal Visual System]. The white and grey circles on each node indicate the level of support values in ultrafast bootstrapping and Bayesian posterior probability. The labels at the branch tips indicate the number of collapsed sequences and the region in which they were sampled. The species names in the dashed boxes indicate the original species name before the revision if different from the name after revision. Detailed information on each sample and 95% CI can be found in the Supporting Information (Table S1; Fig. S3).

Figure 3.

Maximum likelihood phylogeny of the genus Perpolita inferred from 1 171 157 bp that was obtained from ddRAD-seq. The species names in the dashed boxes indicate the original species name before the revision if different from the name after revision. Detailed information on each sample and 95% CI can be found in Supporting Information, Table S1.

The mitochondrial Cytb phylogeny reveals the presence of three well-supported monophylies (UFB > 95/BPP = 1.00; P. binneyana, P. electrina, and P. petronella), and one partially supported monophyly (96/0.52; P. hammonis; Fig. 2; Supporting Information, Fig. S3). Perpolita binneyana occidentalis did not exist as even a poorly defined subclade within P. binneyana (Fig. 2; Supporting Information, Fig. S3). The lower support for P. hammonis is at least partially due to the presence of two intra-specific subclades, with subclade A being well supported (96/1.00) and ranging across Eurasia west of the Urals, and with subclade B being partially supported (98/0.72) and restricted to the Alps except for one haplotype from the Czech Republic (Fig. 2; Supporting Information, Fig. S3).

The largest discrepancy between prior taxonomic concepts and mtDNA sequence was observed in the highly supported species-level monophyly (100/1.00) that had been previously subsumed either into P. hammonis in Eurasia or P. electrina in Alaska. Based on the nomenclatural record—and a cursory examination of potential types—we suspect that Perpolita radiatella (Reinhardt, 1877) is the earliest available nomen for this entity, and we will use it to demarcate this clade throughout the rest of the paper. This entity is again divided into two subgroups, with subclade A (95/1.00) ranging from central and northern Europe to central Siberia, and subclade B (75/0.76) ranging from central Siberia to Alaska (Fig. 2; Supporting Information, Fig. S2). These subclades co-occur within the Altai Mountains of south-central Siberia.

The concatenated ITS1+ELAV8 phylogeny basically replicates the mtDNA topology (Fig. 2; Supporting Information, Fig. S3) with P. dalliana representing a very long branch separated from the remaining Holarctic species. Perpolita binneyana, P. electrina, P. petronella, and P. radiatella formed maximally-supported monophylies (100/1.0). Again, subclade structure in P. binneyana was not geographically correlated, while P. radiatella possessed a strongly supported (99/1.00) geographically restricted subclade ranging west of the Urals. Perpolita hammonis again demonstrated considerable within-clade diversity, with a highly-supported (100/1.0) Alpine subclade again being present.

The ddRAD phylogeny demonstrated a very similar topology but with much greater node support and resolution. Again, P. dalliana existed as a very long branch greatly removed from the remaining species. Within the Holarctic species, P. binneyana, P. electrina, P. hammonis, and P. radiatella each represented a distinct, maximally-supported (100/1.0) monophyly. While the use of only a single individual of P. petronella prevented species-scale clade resolution, its branch length was similar to other species. Perpolita hammonis was found to exist as two fully-supported (100/1.00) subclades, one ranging across Europe and the other being restricted to the Alps. Perpolita radiatella from central Siberia to Alaska again was found to exist in a separate subclade (100/1.00) from European populations. Again, the basal node for Holarctic taxa demarcated P. radiatella from the remaining species.

Divergence time

RelTime analysis (Fig. 4) suggested divergence times ranging from 18.7 Mya for P. radiatella (note that this is also the calibration node) to 13.4 Mya for P. electrina, 11.0 Mya for P. binneyana, and 10.0 Mya for P. petronella and P. hammonis. The divergence time between eastern and western subclades of P. radiatella was estimated at 4.4 Mya and between Alpine and pan-European P. hammonis subclades at 4.0 Mya. BPP analysis generated similar species divergence time estimates of 17.5, 13.4, 10.2, and 10.0 Mya, respectively.

Figure 4.

Divergence time estimation of the genus Perpolita on the maximum likelihood tree. The red circle indicates the node calibrated by the oldest known Perpolita fossil, and the integer number by the node referred to the node number. The blue bars at each node indicate 95% CI estimated with RelTime. The species names in the dashed boxes indicate the original species name before the revision if different from the name after revision. Detailed information on each sample and 95% CI can be found in Supporting Information (Tables S1, S3).

Shell morphology

We found diagnostic shell differences to exist between all genetically-supported Holarctic Perpolita species (Table 4; Fig. 5). The development of sillons (parallel spiral grooves cut into the shell) on the lower surface (Fig. 6; Supporting Information, Fig. S4), especially near the umbilicus, represents a previously unrecognized diagnostic trait. Perpolita binneyana is distinguished by its pale shell colour, low spire, wide umbilicus, and absent sillons. Perpolita petronella also has a pale shell, but it possesses a high spire, narrow umbilicus, lower whorl expansion rate, and weak sillons. Perpolita hammonis has the darkest shell colour of all species, as well as a wide umbilicus, more rapid whorl expansion, and strong sillons. The neotype for this taxa (Oslo Museum D27826) was found to possess all of these features, except paler colour via shell aging. Perpolita electrina has a lighter brown shell colour with a narrower umbilicus, a larger shell with a lower whorl expansion rate, and missing/weak sillons. It also has the widest apertural opening of all species. Perpolita radiatella has a light brown shell with a low spire, moderately wide umbilicus, rapid whorl expansion rate, and weak/missing sillons. It can be distinguished from P. hammonis by much weaker sillons and somewhat lighter shell colour, and from P. electrina by its more open umbilicus, smaller median shell size, and greater whorl expansion rate.

Figure 5.

Shell images of each Perpolita species from genetically validated populations (see Supporting Information, Table S1 for location details defined by sample code): (a) P. binneyana (H216_1; Alaska), (b) P. electrina (H207_1; Main), (c) P. radiatella (H284_1; Japan), (d) P. hammonis (H382_1; France), (e) P. petronella (H040_1; Altai Republic).

Figure 6.

Development of shell bottom sillons: (a–c) P. radiatella; (a) H390, Belarus; (b) H340, Alaska; (c) H346, Czechia; (d–f) P. hammonis; (d) H046, Switzerland; (e) H382, France; (f) H432, Latvia; (g) P. petronella, H038, Norway; (h) P. electrina, H207, Ontario; (i) P. binneyana, H216, Alaska. See Supporting Information, Figure S4 for detail of P. radiatella and P. hammonis shell bottom microsculpture.

Table 4.

Perpolita species shell characters (following molecular revision). See also Figures 5–6 and Supporting Information, Figure S4

Species	Median width and range (mm)	Median height and range (mm)	Colour	Sillons	Spire	Aperture	Umbilicus	Maximum numbers of whorls	Whorl expansion rate (median)
P. binneyana	3.3 (2.8–3.9)	1.7 (1.4–2.0)	Pale ivory	Missing	Low-high	Open	Widely open	3.50	2.0
P. electrina	3.6 (3.2–4.2)	2.0 (1.7–2.2)	Light brown	Missing-weak	Low-high	Widely open	Narrow	3.35	1.9
P. hammonis	3.3 (2.9–4.4)	1.7 (1.4–2.2)	Dark brown	Strong	Low	Open	Widely open	3.55	2.2
P. petronella	3.9 (3.2–4.7)	2.1 (1.7–2.7)	Pale ivory	Weak	High	Open	Narrow	3.70	1.8
P. radiatella	3.3 (2.8–4.1)	1.7 (1.5–2.2)	Light brown	Missing-weak	Low	Open	Open	3.35	2.2

Geographic and ecological range

Each validated species also exhibits a unique geographical and ecological range (Fig. 7). Perpolita hammonis extends from the North Atlantic east to the Urals, and tends to favour moderately humid and neutral-acidic forests but also occurs across a range of wooded and open wetland habitats (Welter-Schultes 2012). Perpolita petronella extends from the Atlantic to Altai Mountains, and appears to favour cooler and more humid sites than P. hammonis, particularly in Europe. Perpolita binneyana extends across the North American taiga from Newfoundland and Labrador to Alaska, south to the Ohio River and down the Rocky Mountains to Arizona, New Mexico, and Texas. It tends to favour upland forest sites. Perpolita electrina ranges in North America from the central and north Atlantic west to the Pacific and north into the taiga, although it appears to be absent from the Beringian parts of Alaska and the Yukon. It is much more of a wetland species than P. binneyana, and is characteristic of fens, marshes, and nutrient-rich wet forests. Both it and P. binneyana are commonly sympatric in wet-mesic taiga. Perpolita radiatella ranges from the Atlantic coast of Scandinavia and central Europe east to Alaska. While sharing a similar environmental niche with P. hammonis, towards the west it tends to be limited to moister sites while it becomes more of an upland forest species in the far east.

Figure 7.

Potential distribution probability of each Perpolita species estimated by ecological niche modelling using 1.3-9 in R (Hijmans et al. 2021). The best models for each species, selected by the area under the curve based on cross-validation were used for the estimates. Occurrences on the map are environmentally pruned samples for modelling. The estimated ice sheets at the LGM are based on Ehlers et al. (2011) and the LGM climate is based on CCSM4.

Habitats supporting sympatric Perpolita (Supporting Information, Table S2) occur throughout the Holarctic and include Scandinavian and eastern European mires and humid forests (P. radiatella/P. hammonis—nine sites, P. radiatella/P. petronella—61, P. hammonis/P. petronella—21, P. radiatella/P. hammonis/P. petronella—nine), central Asian hemi-boreal taiga (P. radiatella/P. petronella—52 sites), and North American cool-temperate forests and taiga from Labrador to Alaska south to Maine, Iowa, and British Columbia (P. binneyana/P. electrina—186; P. binneyana/P. radiatella—two).

Ecological niche modelling

The best Maxent climatic models based on AUC were also the highest-ranking models based on AICc (Table 3; Supporting Information, Table S4). Across all five species and all three climatic circulation models, the estimated appropriate climate area during the LGM was smaller than in the modern era (Fig. 7; Supporting Information, Figs S5, S6). While Perpolita hammonis was isolated to discrete areas northwest and southeast of the Alps and to the northwest coast of continental Europe; P. petronella had a large and more uniform suitable climate region centred in central Europe between the Fennoscandian/British and Alpine ice sheets; P. binneyana and P. electrina existed as overlapping uniform potential distributions both ranging across southeastern North America from the Ozark Mountains to the Carolinas; P. radiatella occurred in approximately three isolated appropriate climate areas, one centred on Europe, one east of the Tibetan Plateau, and the other in southern Alaska. While an appropriate climate for P. radiatella also occurred south of the North American ice sheet, based on modern phylogeography and fossil records it appears this region was never occupied (Fig. 7; Supporting Information, Figs S5, S6).

Discussion

Species diversity in Perpolita

Because the North American Gulf Coast P. dalliana and P. suzannae were found to occupy a distinct very long branch across mtDNA, nDNA, and SNP data, it seems likely they reside within a different sister genus. Because our data shows that these two taxa possess almost identical DNA sequences it seems likely that they are best considered a single species. Based on overall shell appearance and biogeographical range, ‘Helix’ subhyalina L. Pfeiffer, 1867 of southeastern Mexico is likely to be part of the same group with ‘P’. dalliana. ‘Glyphyalus’ quillensis de Winter, van Leeuwen & Hovestadt, 2016, from the Dutch Caribbean island of St. Eustatius, possesses not only a very similar shell, but 95% similarity in the CO1 sequence to P. daliana (De Winter et al. 2016). Given that true Glyphyalina (Glyphyalus) from the southeast USA is genetically very distinct (J. Slapcinsky, personal communication) it is clear that ‘G’. quillensis is part of the same genus-level group as ‘P’. dalliana. This putative new sister genus to Perpolita thus appears likely to range across the central and southern Gulf coast to at least the eastern Caribbean.

Our phylogenetic analyses of the remainder of Perpolita within the temperate/boreal Holarctic consistently show existence of five distinct species-level clades: P. binneyana, P. electrina, P. hammonis, P. petronella, and P. radiatella. Because each possesses a unique suite of conchological traits, geographic ranges, and ecological extents (Table 4; Figs 5–6; Supporting Information, Fig. S4) they can be accurately identified in the absence of molecular data. Only P. petronella fully followed traditional taxonomic concepts and diagnostic features (e.g. Welter-Schultes 2012, Horsák et al. 2013). In the remaining taxa significant errors existed between traditional concepts and empirical genetic data. While P. binneyana was confirmed as a valid species, we found no evidence of gene pool partitioning between eastern and western North American populations. We thus recommend that subspecies distinctions (e.g. P. b. occidentalis) be avoided. Perpolita electrina was found to possess a narrower range than previously thought, being limited from the Atlantic Ocean to British Columbia. Material west from the Yukon actually represented P. radiatella. Perpolita hammonis was found to possess a much more restricted range than traditionally understood, being limited from the North Atlantic (west to Iceland and the Azores) to eastern Europe. All samples initially identified as P. hammonis from the southern Urals eastward, as well as considerable material extending west to Scandinavia and central Europe, actually represents P. radiatella.

As has been previously seen (Nekola and Horsák 2022), errors in traditional taxonomic concepts represent a mix of over splitting, over lumping, and incorrect designation of diagnostic features. These can lead to significant bias in apparent ecological and biogeographic pattern. The existence of these issues within Perpolita is a reminder that traditional taxonomic concepts must be considered untested hypotheses requiring empirical confrontation before being accepted and used in ecological, biogeographic, and biodiversity analyses (e.g. Nekola et al. 2015, 2018, Aksenova et al. 2017, Horsáková et al. 2020, Bespalaya et al. 2023).

What makes this case remarkable is that P. radiatella was overlumped into not one but two separate taxa. The ultimate cause for this is related to the distribution of 18^th and 19^th century taxonomic schools whose ranges were formed when workers had limited access to information from distant areas (Graf 2007, Vinarski 2018). Thus, members of the European School overlumped Eurasian populations of P. radiatella into the European P. hammonis while those in the North American School overlumped Alaska/Yukon P. radiatella populations into North American P. electrina. This led to the generation of a perceived and artificially abrupt biogeographic transition between the North American and Eurasian Perpolita faunas centred on the Bering Strait. In fact, this transition is gradual and spread across the whole of Beringia. Similar spatially-autocorrelated taxonomic error has also been observed in freshwater mollusc faunas between the Russian and Western European schools (Korniushin 1998, Graf 2007, Vinarski and Kramarenko 2015, Vinarski 2018).

Unlike many other land snail genera (e.g. Ainohelix and Ezohelix: Morii et al. 2015, Bradybaena: Hirano et al. 2019, Euconulus: Horsáková et al. 2019, Pupilla: Nekola et al. 2015, Haase et al. 2021; Schileykula: Harl et al. 2020) we observed no topological incongruence in individual clade membership between mtDNA and nDNA data streams. SNP-based divergence times between the five Perpolita species were estimated to be more than 10 Mya. While this is older than many other land snail genera (e.g. Köhler and Criscione 2014, Hirano et al. 2019, Brozzo et al. 2020, Hwang et al. 2021, Neiber et al. 2021), most of these other taxa are limited to isolated oceanic or edaphic islands. The few studies conducted on land snails possessing continental distributions demonstrate divergence times closer to our estimates (Hausdorf and Neiber 2022, Neiber et al. 2022). Thus, our evidence suggests that even though Perpolita species reside within northern areas greatly impacted by cyclic Pleistocene glaciations (Murray-Wallace 2023), this dynamism is likely not responsible for species-scale diversification and cladogenesis. Additionally, distribution of similar shell colours and habitat requirements across the phylograms suggests that intercontinental species pairs sharing these traits may have independently evolved these features, and that they do not represent ancestral traits.

Intraspecific diversity in Perpolita

Considerable variation in intraspecific gene pool size was noted between the five Perpolita species, ranging from little in P. binneyana, P. electrina, and P. petronella, to two principal subclades (likely representing separate subspecies) in P. radiatella, to a complex multi-clade structure in P. hammonis. Modern range size was a poor predictor with smallest modern range taxa (P. hammonis) possessing the most complex gene pool.

However, climate modelling of potential LGM distributions appears to explain these levels. Species which were predicted to possess large continuous LGM ranges (P. binneyana, P. electrina, P. petronella) lack clear geographic structuring in their modern gene pools. Perpolita radiatella possesses two well-supported intraspecific clades, one currently ranging from central Europe to central Asia, and the other from central Asia to Alaska. Both coexist within the Altai Mountains of south-central Siberia. Projected LGM distribution for P. radiatella suggests a highly disjunct distribution, with distributional centres located in central-east Europe and eastern Beringia. It is interesting to note that similar European vs. Beringian subpopulations have been reported (Nekola et al. 2015, Horsáková et al. 2018) in the land snails Pupilla alpicola and Euconulus alderi, with the latter harbouring disjunct European vs. Beringian projected LGM ranges (Nekola et al. unpubl. data).

In spite of its possession of the most limited modern range of all Perpolita species, P. hammonis demonstrated the most complex intraspecific gene structure: not only are there two well-supported subclades (one ranging across the entire species range and the other being limited to the Alps) but there are a number of additional more poorly-supported apparently spatially-restricted groups as well. LGM range reconstructions suggest its occurrences were fragmented into multiple centres of occurrence ranging from the North Atlantic coast to the west, south, and east of the Alps. The positive relationship between LGM refugia number and modern gene pool size has been previously noted in other taxa groups (e.g. Cheddadi et al. 2006, Dussex et al. 2014, Roberts and Hamann 2015). These results suggest that LGM climate change may underlie levels of subspecific diversification in Perpolita.

Acknowledgements

We are grateful to Stefan Meng for providing the materials of Perpolita from eastern Siberia. We also thank Yasuto Ishii, Satoshi Chiba, Shu Nishida, and Isao Sano for supporting the ddRAD work, and Jan Divíšek for advice on ENM methods.

Author contributions

M.H., J.C.N., and V.H. collected the samples; M.H. and J.C.N. conceive the idea; T.S., M.N., E.L., and T.H. conducted the molecular work with partial support of V.H.; T.S. analysed the data and led the writing with the support of J.C.N. and M.H.; all authors reviewed and approved the final version of the manuscript.

Conflict of interest

None declared.

Funding

The study was financially supported by the Czech Science Foundation (22-230055) and the dd-RAD work was financially supported by the Japan Society for the Promotion of Science (JSPS), Grants-in-Aid for Scientific Research (KAKENHI: 23K14261), and the Showa Seitoku Memorial Foundation (Reiwa-4).

Ethics and Permits

All research pertaining to this article did not require any research permits.

Data availability

All data generated during this study are included in this published article, its supplementary information files and deposited in GenBank (PP565243–PP565334, PP573309–PP573503 & SRR28509296-SRR28509312).