The primary material for this study was collected between 5 and 21 July 2014 by two of the authors (M.V. and I.N.) in several localities situated on the western coast of Greenland. The sampled habitats were concentrated in three areas (Fig. 1):

• • vicinities of Kangerlussuak Settlement (approximately 67°00′N; 50°41′W). The lymnaeids were sampled from 11 waterbodies belonging to three different types (see below);
• • Nuuk Area. We inspected six waterbodies situated in Nuuk (the capital of Greenland) or in mountain ridges surrounding the town. The lymnaeids were found in a single waterbody – a small shallow lake lying in the center of Nuuk (64°11′15′′N; 51°43′14′′W);
• • vicinities of Narsarsuaq settlement. More than 10 waterbodies of various types were inspected. Lymnaeids were found only in few lakes situated near Quarsiarsuq settlement (61°08′09′′N; 45°32′17′′W). In other waterbodies snails were represented by the planorbid species Gyraulus arcticus (Møller, 1842) only, or were absent;

In the vicinities of Kangerlussuak, three types of waterbodies could be discerned (Fig. 2):

• • large mountain lakes with heavy wave action and stony substrate. Lake Tasersuatsiaq (Ferguson) is a typical representative of this group (see Fig. 2A);
• • small tundra lakes with soft bottom substrate and a narrow belt of vegetation. This is the most abundant type of lake in the area, most of them are nameless (see Fig. 2B);
• • floodplain pools. Very shallow, with no vegetation and presumably non-permanent (see Fig. 2C). We inspected several such pools situated in the floodplain of Qinnguata Kuussua River outflowing from the Russell Glacier.

Snails were picked up by hand from the substrate in the shallow zones of lakes or from vegetation. Also, a sieve was used to collect mollusks from the bottom. The animals sampled were immediately put into plastic tubes filled with 96% ethanol. In total, more than 1500 specimens were collected (most of them were empty shells abundantly sampled for conchometric study). The bulk of the sample is kept in the collection of the Laboratory of Macroecology and Biogeography of Invertebrates, Saint-Petersburg State University, Russia. The shells and tissues of individuals used for genetic studies are in the collection of Senckenberg Natural History Collections Dresden (Museum of Zoology), Germany (SNSD hereafter). Some voucher specimens were given to the collections of the Zoological Institute of the Russian Academy of Sciences, Saint-Petersburg, and Vienna Museum of Natural History, Austria (NHMW hereafter).

The type series of three nominal Møller's [19] species of Greenland Lymnaeidae are kept in the Natural History Museum of Denmark (Copenhagen, ZMUC hereafter) and were examined by us prior to the field work in West Greenland. Several samples of Greenland lymnaeids kept in NHMW were also examined.

The type specimens of Limnaea vahlii Møller, 1842 are kept in ZMUC under the numbers GAS-106, 107, 108, 109, and 110. There are several hundreds of syntypes in ZMUC and six others are housed in the Stockholm Museum of Natural History, Sweden (No. SMNH 3842). One of the syntypes of L. vahlii was previously illustrated by Schiøtte and Warén [18, their fig. 7]. The type locality of L. vahlii is not known exactly. According to Schiøtte and Warén [18], the syntypes were collected in ponds of southwestern Greenland, without exact localisation of these waterbodies.

The type series of Limnaea holboelli Møller, 1842 is kept in ZMUC under the numbers GAS-54, 55, and 56, and includes 27 syntypes collected somewhere in West Greenland. One of the syntypes was illustrated by Schiøtte and Warén [18, their fig. 8].

The type series of Limnaea pingeli Møller, 1842 (ZMUC No. GAS-83) includes five syntypes. The type locality of this species is “a small pond between the mountains at Godthaab” [18, p. 12]. Godthaab is the former name for Nuuk (see above).

Shells of lymnaeids were measured following the standard scheme [11,26] by means of calipers of the ocular-micrometer of a stereoscopic microscope with an accuracy of 0.1 mm. Six measurements were taken from each shell: shell height and width, spire height, body whorl height, height and width of the aperture. The number of whorls of each shell was also counted. The measurements were then analyzed statistically by means of the multivariate statistical algorithms. Discriminant analysis and analysis of principal components (PCA) were realized with the STATISTICA 6.0 software (StatSoft Inc., USA).

The application of different multivariate techniques to the studies of variation in shells of freshwater gastropods is usually complicated by inevitable size differences among specimens in large samples. There are several ways how to set aside this factor and to analyze the “pure” diversity in shape [27]. In this study, we used PCA to extract several principal components, the first of which (PC1) is typically tightly correlated with absolute measurements of shells and thus is considered as the composite “size” variable [28,29]. Then, we carried out the discriminant analysis based on the values of PCs 2–4, which are presumably less correlated with absolute dimensions and reflect mainly differences in shape [29].

Dissections of snails were made following the standard procedure [11]. We examined the taxonomically significant parts of the male and female reproductive systems, including the copulatory apparatus, prostate, bursa copulatrix, and its duct. The index of the copulatory apparatus (ICA hereafter) was calculated for each dissected individual as the ratio between of the lengths of its praeputium and penis sheath (see Fig. 7C and [30] for details).

In total, 176 shells (including syntypes) were measured and 64 snails were dissected (Table 1). Unfortunately, many of the type specimens of Møller's taxa are damaged (broken shell apices) and thus were unavailable for measurement. The juvenile and subadult shells were excluded from morphometric analysis as well. Therefore, the overall number of type specimens used in this study is relatively low compared to the total number (see Table 1).

For outgroup comparison in the molecular genetic analyses we included sequences of Palearctic specimens of the species Aplexa hypnorum (Linnaeus, 1758) from the freshwater gastropod family Physidae and from the lymnaeid species Galba truncatula (O.F. Müller, 1774), Lymnaea stagnalis (Linnaeus, 1758), L. taurica (Clessin, 1880), Stagnicola corvus (Gmelin, 1791), S. palustris (O.F. Müller, 1774) and S. fuscus (C. Pfeiffer, 1821). We included own sequences of Ladislavella terebra (Westerlund, 1885), L. elodes (Say, 1821), and L. catascopium in the in-group. In addition, we used ITS2 sequences of L. elodes, L. catascopium, L. terebra [= Catascopia occulta (Jackiewicz, 1959)] and L. emarginata (Say, 1821) available from GenBank.

Tissue samples taken from the foot were fixed in 100% ethanol or 100% isopropyl alcohol. They were registered in the tissue collection of the SNSD with both a new tissue voucher number and the collection number in the mollusk collection of SNSD and stored at –80 °C.

For the molecular genetic analyses, we obtained sequence data of the nuclear ITS2 spacer (partial, with 28S ribosomal RNA gene, partial sequence, 280 bp in A. hypnorum, up to 506 bp in L. taurica, the length of the ITS2 spacer varies within the genera and families) and a 329 bp fragment of the cyt-b gene as mitochondrial marker. For primers and protocols of DNA extraction, Polymerase Chain Reaction (PCR), purification of PCR products and DNA sequencing, see Schniebs et al. [31].

All specimens used in the molecular genetic study are listed in Appendix 1. All DNA-sequences have been placed in the European Nucleotide Archive (ENA, see http://www.ebi.ac.uk/ena/).

The alignment of the cyt-b gene dataset was performed using the ClustalW algorithm implemented in MEGA6 [32]. ITS2 sequences were also aligned using the ClustalW algorithm of MEGA6 [32], because of the lack of any analogues with known secondary structure in the ITS2 database [33]. To eliminate highly variable fragments from the ITS2 alignment, we used GBlocks v. 0.91b [34,35] using options for less stringent selection enabling inclusion of gap positions, smaller final blocks and less strict flanking positions. The resulting length of the ITS2 alignment was 468 bp (85% of the original 545 positions). Each alignment was collapsed into unique haplotypes (see Appendix 1) using an online FASTA sequence toolbox (FaBox 1.41; [36]). Absent sites were treated as missing data. A partition homogeneity test was applied in PAUP v. 4b10 to confirm the congruence of phylogenetic signals among sequence data sets [37,38]. This test reveals significant conflict (P = 0.01), and we therefore analyzed each gene data set separately. Each codon of the cyt-b gene was treated as a separate partition.

All phylogenetic analyses were performed at the San Diego Supercomputer Center via the CIPRES Science Gateway [39]. Maximum likelihood (ML) analyses were conducted using RAxML v. 8.2.6 HPC Black Box [40]. A unique GTR model was applied for each partition with corrections for gamma distribution. Nodal support values were estimated using an automatic, rapid bootstrapping algorithm according to the developer's recommendation [41], and the majority-rule consensus tree was constructed from the independent searches. Bayesian inference (BI) analyses were performed in MrBayes v. 3.2.6 [42]. The best models of sequence evolution for each partition based on the corrected Akaike Information Criterion (AICc) of MEGA6 [32] are presented in Table 2. Two runs, each with three heated (temperature = 0.1) and one cold Markov chain, were conducted for 30 million generations. Trees were sampled every 1000th generation. After completion of the Markov Chain Monte Carlo (MCMC) analysis, the first 15% of trees were discarded as burn-in (pre-convergence part), and the majority-rule consensus tree was calculated from the remaining trees. Convergence of the MCMC chains to a stationary distribution was checked visually based on the plotted posterior estimates using an MCMC trace analysis tool (Tracer v. 1.6; [43]). The effective sample size (ESS) value for each parameter was always recorded as > 10,000.

To delimit putative species-level units, we used a molecular approach based on the concept of Molecular Operational Taxonomic Units (MOTUs thereafter) [44,45]. MOTUs were separated using the Poisson Tree Processes (PTP) model to infer putative species boundaries on a phylogenetic input tree inferred from a ML analysis of the cyt-b haplotype sequences [46]. We used a Bayesian implementation of the PTP model for species delimitation through an online bPTP server (http://species.h-its.org/ptp) with 500,000 MCMC generations and 15% burn-in [46]. The outgroup haplotypes were removed from the input tree using an appropriate option of the server.

To test the null hypothesis of equal evolutionary rate throughout the tree (a global molecular clock model), we performed the likelihood ratio test (LRT) [32,47] based on the consensus phylogenetic tree inferred from ML analyses of the cyt-b haplotypes. The null hypothesis was rejected at a 5% significance level (P = 0). Moreover, the external molecular clock rates appear to be hardly applicable to the Lymnaeidae, because of the possibility of ultrafast molecular evolution in several taxa [9]. The fossil-calibrated phylogenetic reconstructions for lymnaeid snails are still lacking. For biogeographic reconstructions, we therefore used an uncalibrated ultrametric tree of the cyt-b haplotype sequences, which was obtained using BEAST v. 1.8.4 with a lognormal relaxed clock algorithm and the Yule speciation process as the tree prior [48]. Calculations were performed at the San Diego Supercomputer Center through the CIPRES Science Gateway [39]. We specified the same substitution models to the three-codon positions as in the MrBayes analyses, but have used the simple HKY model instead GTR. Two replicate searches were conducted, each with 30 million generations. The log files were checked visually with Tracer v. 1.6 for an assessment of the convergence of the MCMC chains and the effective sample size of the parameters [43]. Almost all ESS values were recorded > 1000; the posterior distributions were similar to the prior distributions. The resulting tree files from two independent analyses were compiled with LogCombiner v. 1.8.4 [48]. The first 15% of trees were discarded as an appropriate burn-in (pre-convergence part). The maximum clade credibility tree was obtained using TreeAnnotator v. 1.8.4 [48]. We tested ancestral area patterns using three different approaches, i.e., Statistical Dispersal-Vicariance Analysis (S-DIVA), Dispersal-Extinction Cladogenesis (Lagrange configurator, DEC), and Statistical Dispersal-Extinction Cladogenesis (a nonparametric ‘Bayes–Lagrange’ approach to biogeographic inference, S-DEC) implemented in RASP v. 3.2 [49]. For an ancestral area reconstruction, we used the set of 18,002 uncalibrated trees that were combined from two runs of BEAST v. 1.8.4 with an appropriate burn-in (see above). As a condensed tree, we used the user-specified consensus tree, which was obtained based on this set of trees using TreeAnnotator v. 1.8.4 (see above). Only Ladislavella taxa were used in the analyses, all the other taxa were removed from input trees using the appropriate option of RASP v. 3.2. We coded three distribution areas as follows: Greenland (A), mainland North America (B), and Eurasia (C). The S-DIVA models were calculated with the following parameters: max areas = 3; allow reconstruction with max reconstructions = 100 and max reconstructions for final tree = 1000; allow extinctions. The DEC and S-DEC analyses were run with default settings and max areas = 3. In addition to the evaluations obtained from each analysis separately, we used the generalized results of all three modeling approaches that were combined using an algorithm implemented in RASP v. 3.2.

The variability of shell size and proportions among syntypes of Lymnaea vahlii proved to be very high. The shell shape varies from tall and slender turriculate shells to almost ovate specimens with short conical spires (Fig. 3A–D). Most syntypes are brown coloured, fragile, with thin shell walls, and with more or less prominent axial lines on the body whorl surface (see Fig. 3C and D). The aperture proportions are also very variable. Despite the obvious difference among the shells of syntypes, all of them more or less correspond to Møller's [19] original diagnosis, which, however, is very brief and vague.

The extent of conchological variation in the L. holboelli syntypes is comparable with that in the type series of L. vahlii. The shells are diverse in their size and especially in their shape (see Fig. 3E–G). It proved impossible to discriminate the two species based on their shell proportions. Møller [19] believed that L. holboelli shell has a more acute spire and a deeper suture compared to L. vahlii. We should note, though, that these differences appear not to be stable, especially if the overall extent of conchological variation within the type series of L. vahlii is considered. The conchological variables of the two species greatly overlap (Table 3). The two shells of the syntypes of L. pingelii are rather similar to the syntypes of the two above-mentioned species (see [18], their Fig. 9) and, most probably, rather represent an intraspecific form of L. vahlii.

According to our field observations and shell measurements, the intrapopulation variation in the Greenland lymnaeids is very high, and its amplitude corresponds to the extent of shell variation found within the type series of L. vahlii. Nevertheless, one may note that the shells can be sorted into two more or less prominent groups differing in their size and proportions. The first group includes snails with tall and slender high-spired shells, resembling shells of species belonging to the Stagnicola elodes group sensu Burch [50], which are widely distributed in North America (compare Fig. 4A–C and D–E). This group proved to be 100% distinguishable by means of the discriminant analysis (Table 4). The shells of the second group are ovate-conical or almost ovoid, with short spires and strongly inflated body whorls. They look similar to shells of species of the Stagnicola emarginata/catascopium group sensu Burch [50], also of Nearctic origin (compare Fig. 4F–H and I–K). The two forms will be designated hereafter as the ‘narrow’ and ‘broad’, respectively. These can be characterized mainly by differences in two conchological ratios: shell width:shell height and body whorl height:shell height (see Table 2).

In the vicinities of Kangerlussuaq, snails of the ‘narrow’ form occur almost exclusively in small tundra lakes and floodplain pools, whereas specimens belonging to the ‘broad’ group were found in two large mountain lakes in vicinities of Kangerlussuak settlement: Tasersuatsiaq (Ferguson) and Saningassoq (only eight specimens were collected here). The possible topotypes of L. pingelii sampled in the city of Nuuk also belonged to the ‘broad’ form (see Fig. 4H), although their habitat was a small lake. However, this waterbody lies in a small depression between two mountain gorges, and its bottom was covered by large and little stones (see Fig. 2D). Thus, it is likely that the distribution of both morphotypes is not random, but rather depends on the type of waterbody, but our evidence is insufficient to establish this pattern with certainty.

We found a difference between ‘broad’ and ‘narrow’ forms in the absolute size of their shells. The ‘broad’ shells are smaller and this size difference is statistically significant (Student's t-test: t = 5.864; P = 0.00). Therefore, the direct comparison of the two forms is difficult, unless the size factor of variation is removed. This was achieved with the PCA analysis of conchological variables (see Material and methods). Examination of the axes’ loadings reveals that the first axis of PC1 (explaining 79.0% of the variation) is the ‘size’ variable, since all linear measurements tightly correlate with PC1 values (Table5). The values of the subsequent axes (PCs 2–4) correlate with the number of whorls, spire height, and aperture height; hence, these variables contribute most to the differences between ‘narrow’ and ‘broad’ conchological forms.

The discriminant analysis performed after setting aside the first axis still allowed us to distinguish between the ‘broad’ and ‘narrow’ specimens, although the overall exactness of discrimination dropped to 80.0% Table 6. The phenotypic distinctness is not complete due to the existence of ‘intermediate’ specimens. The ‘broad’ and ‘narrow’ individuals form a single ‘cloud’ of points in the multivariate space (Fig. 5) that indicates the absence of the full hiatus between them (see also Table 2). The exactness of the discrimination between ‘narrow’ snails living in floodplain pools and tundra lakes (waterbody types II and III) is rather low (for ‘narrow’ specimens from floodplain pools, it is equal to 57.5%).

The direct comparison between specimens collected in 2014 and the syntypes of Møller's nominal species was not possible due to large differences in shell size: the syntype shells are roughly two times smaller than the newly sampled individuals (see Table 2). Therefore, the procedure of setting aside the size component of variation was repeated on a united sample consisting of syntypes and non-syntypes. The analysis was successful in eliminating the size variation, since no differences between the type specimens and newly collected shells were found. All syntypes clustered together with the ‘narrow’ specimens sampled in 2014 (Fig. 6).

The structure of the copulatory organ in the Greenland lymnaeids is uniform across populations and waterbody types. All snails have an oblong, usually dark-pigmented praeputium with a more or less prominent distal end (near its border to the penis sheath). The penis sheath is much shorter and narrower than the praeputium, and strongly inflated at its distal end, forming a club-like structure (Fig. 7). The ratio between the lengths of the praeputium and the penis sheath (ICA) varies from 1.26 to 2.32. Though the mean ICA values vary across populations (Table 7), the extent of their variation corresponds to those observed in other lymnaeid species [30] and thus this trait is not effective for discrimination of populations and conchologically defined forms of L. vahlii.

The bursa (receptaculum seminis) is short and compact, its duct is rather long and very narrow (see Fig. 7G).

Though we had no large samples of other Nearctic species of Ladislavella to carry out a comparative analysis, the examination of the internal structure of several individuals of L. elodes from USA waterbodies revealed no substantial differences in the structure of the copulatory organs and the female reproductive organs (see Fig. 7).

The mantle pigmentation of the Greenland Ladislavella proved to be variable as well (Fig. 8). Individual snails differ primarily in the size of the black patches and in the pattern of their distribution over the mantle surface. The largest black patch was found in a specimen from Lake Tasersuatsiaq, but, as it will be shown below, this difference does not correlate with genetic distance, making it unsuitable as a basis for any phylogenetic or taxonomical inferences.

The hypotheses for the phylogenetic relationships of Ladislavella vahlii recovered from Bayesian analysis of cyt-b and ITS2 sequences are represented in Fig. 9. Both trees demonstrate a similar topology, with two large groups of lymnaeids found: the first one includes taxa of supposedly Neartic origin (Galba, Ladislavella), whereas the second cluster embraces two Old World genera, Lymnaea and Stagnicola. This dichotomy has full Bayesian posterior probabilities/ML support in both cases.

The results of the species delimitation analysis indicate that the Ladislavella catascopium s. l. clade contains three MOTUs, each of which may represent a separate species-level cryptic taxon (see Fig. 9A). The most divergent haplotypes are from California (46m) and from Illinois (45m), whereas the other haplotypes belong to a single MOTU corresponding to Ladislavella catascopium s. str. (Nos. 40m–44m). The Greenland lymnaeids formed a separate cluster within Ladislavella catascopium clade in both trees. Three mitochondrial haplotypes unique for Greenland were recovered during this study (haplotypes 40m, 41m, and 42m in Fig. 9A). Two specimens of L. catascopium s. str. collected in British Columbia shared the same haplotype (44m), which proved to be the most similar genetically to the Greenland snails. The haplotypes 43m, 45m, and 46m were found in snails from the United States, provisionally determined as L. emarginata (see Appendix 1).

Two unique nuclear haplotypes (42n and 43n) were found in the Greenland lymnaeids. These haplotypes have not been revealed in the mainland populations of L. catascopium s.l. (see Fig. 9B). Seven other nuclear haplotypes (Nos. 39-41 and 44-47, see Fig. 9B) are distributed throughout North America, from British Columbia to California. Their distribution, however, is not spatially patterned: the snails sampled in Michigan do not constitute a separate cluster of the cladogram similar to that formed by the Greenland individuals. The same was characteristic of lymnaeids originated from British Columbia (see Fig. 9B).

The ‘broad’ and ‘narrow’ conchological forms of L. vahlii outlined above proved to be genetically indistinguishable. For example, specimens from Lake Ferguson (belonging to the ‘broad’ form) shared identical haplotypes (40m and 42n) with snails collected from a nameless small lake near Kangerlussuak, which were classified as ‘narrow’ (see Fig. 9, Appendix 1). On the other hand, the distribution of the haplotypes among specimens of the Greenland Ladislavella seems to be governed partly by geographic distances separating sampled areas. The majority of the Greenland haplotypes recovered during this study appeared in a single area and were not found elsewhere. The haplotypes 40m and 42n are characteristic of the Kangerlussuak area, whereas the haplotype 41m was found in Nuuk, and the haplotype 43m – in vicinities of Narsarsuaq settlement. The only exception is the haplotype 43n, whose bearers were collected in all three sampling areas (see Appendix 1). Anyway, we could not find two separate clades that would correspond to the ‘narrow’ and ‘broad’ groups specimens in both our trees.

The existence of a set of haplotypes (both mitochondrial and nuclear) unique for Greenland as well as its isolated position on phylogenetic trees allowed us to consider the Greenland samples as belonging to a separate subspecies indicated in the trees (see Fig. 9A and B) as Ladislavella catascopium vahlii.

It is interesting also to note that all three alleged species of the L. catasopium group living in the North American mainland (i.e. L. catascopium s. str., L. elodes, and L. emarginata) were characterized by their own set of unique haplotypes, not found in congenerics. Thus, L. catascopium s. str. in our study possess the haplotypes 44m, 39n, 44n and 45n; L. emarginata has two own nuclear haplotypes (40n and 41n), whereas L. elodes demonstrates the widest array of haplotypes unique for this taxon, including three mitochondrial (43m, 45m, and 46m) and three nuclear (46n, 47n, 48n) ones.

The combined results of our biogeographic modelling (Fig. 10; Table 8) suggest the Nearctic origin of the L. catascopium clade via a vicariance event separating the former continuous Holarctic range of the L. catascopium–L. terebra MRCA (probability = 0.98). The further radiation of L. catascopium s.l. occurred in North America (probability = 0.95). The origin of the Greenlandic subspecies most likely caused by a relatively recent vicariance event, with subsequent radiation within the island (probability of both events = 1.00).