Phylogeography and systematics of Pyramidula (Pulmonata: Pyramidulidae) in the eastern Alps: still a taxonomic challenge

Abstract

The rock-dwelling land snail Pyramidula pusilla was analysed genetically and morphologically with an emphasis on its eastern Alpine distributional area. Genetic variation and phylogeographic patterns were inferred by mitochondrial cytochrome c oxidase subunit 1 gene sequences (650 bp), while shell variability was investigated by shell measurements and landmark analyses. The phylogenetic analyses revealed two clades in Austria (clades 3 and 4), of which one also contained individuals from Turkey. Other samples originating from the Balkan region (Albania, Greece) and Turkey were well separated in the tree, forming two distinct sister clades (clades 1 and 2). Concerning clades 3 and 4 there is no clear geographical pattern and both clades cooccur at several (11) sites; the same is true for clades 1 and 2. There are also morphological differences between clades: clades 1 and 2 differ from the other clades in size as well as shape, whereas clade 3 differs from clade 4 only in size. Considering all results combined, they imply that the samples analysed might actually represent more than one species and that, contrary to former assumptions, two species of Pyramidula may occur within the eastern Alps. However, comprehensive taxonomic conclusions cannot yet be drawn, for reasons connected with insufficiently justified species assignment: (1) high morphological variation within clades, (2) lack of clear diagnostic information in original descriptions, causing a lack of reliable distributional data and (3) lack of data over the whole distributional range of the genus. Analyses of ncDNA (e.g. microsatellites) to quantify the extent of gene flow, as well as breeding experiments between different clades, could help to clarify the delimitation of species. Despite these taxonomic uncertainties, the results provide insights on Pyramidula within the eastern Alpine region, suggesting that it survived the Pleistocene oscillations without dramatic bottlenecks and is capable of fast expansion to (re)populate formerly glaciated areas.

INTRODUCTION

The rock-dwelling land snail Pyramidula pusilla (Vallot, 1801) belongs to the family Pyramidulidae, of which the sole genus Pyramidula Fitzinger, 1833 is distributed throughout the temperate climatic zone of the Palaearctic from Spain to Japan (Pilsbry & Hirase, 1902). Five of the six European species that are at present accepted in the literature share conchological features, such as a dark yellow to brownish, low-conical and very small shell with a simple, very fragile aperture (Gittenberger & Bank, 1996). Pyramidula chorismenostoma (Westerlund & Blanc, 1879) differs in its scalariform shell.

Pyramidula pusilla is an ovoviviparous species and is the most common and by far the most widespread species of the genus in Europe, distributed in the entire Mediterranean area as well as in Central and Western Europe, where it inhabits sunlit limestone rocks at elevations up to 2,600 m above sea level (asl). It feeds on endolithic lichens (Kerney, Cameron & Jungbluth, 1983; Gittenberger & Bank, 1996). Even though it is described as calciphilous, the species can also be found on rocks containing low amounts of calcium carbonate (Klemm, 1974). Its shells are described as variable in morphology (Pilsbry, 1927–1935), ranging from 2 to 3 mm in size and being broader than high. The four to five finely ridged whorls are connected by a deep suture. While the shell is dark brown in juveniles and young individuals, that of older individuals fades to greyish white (Klemm, 1974; Kerney et al., 1983; Gittenberger & Bank, 1996).

According to Gittenberger & Bank (1996), the distributional area of P. pusilla overlaps with those of nearly all other European species of the genus and, in many parts of Europe, some species even occur syntopically. For example, in the UK P. pusilla cooccurs with P. umbilicata (Montagu, 1803), in Greece with P. chorismenostoma and in northern Dalmatia, western Turkey and Greece with P. cephalonica (Westerlund, 1898) (Gittenberger & Bank, 1996). On the Iberian Peninsula it may cooccur with P. jaenensis (Clessin, 1882) and P. rupestris (Draparnaud, 1801). Moreover, P. pusilla also occurs syntopically with P. rupestris throughout the Mediterranean area to Israel in the southeast (Gittenberger & Bank, 1996).

The taxonomy of Pyramidula changed frequently following the first descriptions of P. pusilla and P. rupestris. Vallot (1801) described P. pusilla as Helix pusilla with unspecified type locality in France. In the same year, Draparnaud (1801) described P. rupestris as H. rupestris, also from ‘France’. After Fitzinger (1833) introduced the genus Pyramidula, P. pusilla was synonymized with P. rupestris by Drouet (1867) until the name was reestablished by Gittenberger & Bank (1996). Generally, four of the six European species (P. pusilla, P. umbilicata, P. cephalonica and P. rupestris) were considered as ‘forms’ or subspecies of P. rupestris before the publication by Gittenberger & Bank (1996). Therefore, one can assume that many older references of P. rupestris in the literature, especially from localities outside the presently known distributional area of P. rupestris, but inside the range of P. pusilla (as described by Gittenberger & Bank 1996), actually apply to P. pusilla. The differences of shell shape and size of all species—except P. cephalonica—are small, so that identifications are usually problematic. The few data in the literature suggest that morphological measurements overlap considerably among species (Gittenberger & Bank, 1996; http://www.animalbase.org/ accessed in November 2014). So far, neither extensive morphometric analyses have been performed on species of Pyramidula, nor have genetic markers been used to provide insights into their phylogenetic relationships. Therefore, the degree of intraspecific variation—genetic as well as morphological—is still unknown. Considerable intraspecific variation and complex phylogeographic patterns are common among landsnail taxa of the northeastern and southeastern Alps, e.g. Arianta arbustorum (Haase et al., 2003), Trochulus hispidus (Duda et al., 2010; Kruckenhauser et al., 2014), Orcula (Harl et al., 2014) and Clausilia dubia (Jaksch, 2012). The Alps feature a wide range of habitat types that differ in altitude, geology and climate. Moreover, Europe was strongly affected by the Pleistocene climatic oscillations, which caused local extinction, isolation in glacial refuges and recolonization events after retreat of the ice sheet. In this regard P. pusilla with its wide distribution range in Europe is of special interest. Because it is quite common in the eastern Alps, especially in formerly glaciated regions (limestone mountain ranges or limestone islands within silicate terrain), the questions arise: where did the recolonization start after the last glacial maximum and do present genetic patterns indicate past glacial bottlenecks?

In the present paper, we attempt a first survey of the genetic variation as well as the shell variability of P. pusilla, with an emphasis on its eastern Alpine distributional area. For the first time the species has been analysed both genetically and morphometrically, to draw conclusions on its phylogeographic history in the eastern Alps. Besides the analysis of shell measurements and landmark data, an c. 650 bp section of the mitochondrial cytochrome c oxidase subunit 1 (COI) gene was sequenced, serving as a basis for tree and network analyses. We asked whether there is any geographic pattern or whether distinct groups can be detected with genetic and/or morphological methods, which might indicate Pleistocene refugia or longterm isolation of populations. This study should provide preliminary data for a comprehensive analysis of the whole genus Pyramidula, leading finally to more clearcut species diagnoses.

MATERIAL AND METHODS

Taxon sampling

The analysed specimens of Pyramidula pusilla were collected from July 2007 to June 2012 from 98 sampling sites. Most of the sampling sites (69 localities: 258 individuals) are located in the Austrian Alps. Specimens from other countries were included as follows: Germany (6 localities: 14 individuals), Turkey (6: 22), Switzerland (5: 13), Italy (3: 9), Slovakia (3: 22), Croatia (2: 8), Slovenia (2: 5), Greece (1: 3) and Albania (1: 3) (Supplementary Material Table S1). The specimens from the last two countries had been identified as P. cf. pusilla. The sampling locations were situated at altitudes between 360 and 2,606 m asl. All specimens collected were preserved in 80% ethanol and stored at 4 °C. For the analyses adult individuals were used preferentially, although the determination of the adult state is somewhat difficult as the aperture of Pyramidula lacks any structures indicating maturity. Each specimen analysed was labelled individually (sample codes consist of an abbreviated taxon name and a four-digit number). All in all, 357 specimens were used for the genetic analyses, 143 of which were also analysed morphometrically (Supplementary Material Table S2). Vouchers of the analysed specimens and their DNA are stored in the DNA and tissue collection of the Natural History Museum Vienna (NHMW).

DNA extraction, amplification and alignment

DNA was extracted from a piece of muscle tissue from the foot with the GEN-IAL First-DNA All-tissue DNA-Kit (GEN-IAL GmbH, Germany), following the manufacturer's protocol. We amplified a section of the mitochondrial COI gene by polymerase chain reaction (PCR) [25 µl reaction volume: 1 µl DNA, 5 µl Q-solution, 2.5 µl buffer, 500 µM dNTPs, 1.5 µl (1.5 mM) Mg²⁺, 0.25 µl of each primer (50 pmol/µl of both COIfolmerfwd 5′-GGTCAACAATCATAAAGATATTGG-3′ and COIschneckrev 5′-TATACTTCTGGATGACCAAAAAATCA-3′; see Duda et al., 2010), 0.1 µl TopTaq-polymerase (0.5 units), 13.9 µl a.d.]. The cycling protocol was: initial denaturation at 94 °C for 3 min, 35 cycles denaturation at 94 °C for 20 s, annealing at 48 °C for 30 s, extension at 72 °C for 60 s and a final extension at 72 °C for 7 min. The PCR products were purified and sequenced in both directions at LGC genomics (Berlin, Germany) and the sequence data was manually aligned and edited with BioEdit v. 7.1.3.0 (Hall, 1999). The alignment used for the analyses had a length of 655 bp. Sequences determined in the course of this study were deposited in GenBank under the accession numbers KT326334 - KT326690.

Phylogenetic analyses

The dataset of 357 sequences of Pyramidula (plus one sequence of Orcula dolium (Draparnaud, 1801) from a sampling site located in the northeastern Alps as outgroup) was used to calculate pairwise p-distances and mean distances between and within the groups/clades with MEGA v. 5 (Koichiro et al., 2011). The program DnaSP v. 5 (Librado & Rozas, 2009) was used to estimate nucleotide and haplotype diversities. We used jModelTest 2 (Darriba et al., 2012) to test for best-fit nucleotide substitution models for tree reconstructions, resulting in selection of the HKY + G + I model. A Bayesian-inference (BI, 5,000,000 generations, default settings) tree was calculated with MrBayes v. 3.2.1 (Huelsenbeck & Ronquist, 2001). Besides the BI tree, we also calculated a Maximum-likelihood (ML) tree using the same dataset with MEGA v. 5 (Koichiro et al., 2011). Network analyses were conducted with the program Network v. 4.6.0.0 (Polzin & Daneschmand, 2003, available at http://www.fluxus-engineering.com) to compute median joining (MJ) networks using the complete COI dataset and also each clade separately, applying the default settings.

Morphological analyses and statistical tests

For the morphological investigation 143 specimens were selected and photographed before they were subjected to the genetic analysis. Pictures were taken of top, bottom and lateral views with a Wild Makroskop M420 and a Nikon DS camera control unit DS-L2, all with the same magnification. We applied both distance-based (traditional) morphometry and shape-variation morphometry based on landmark data (Fig. 1). Since the species cannot be distinguished by the aperture or by any other shell structures, we decided to take the following common measurements: (1) height of each of the first four whorls (H1–H4), (2) width of the protoconch (W0), (3) width of each of the first four whorls (W1–W4), (4) umbilicus diameter (UD) and (5) shell height (SH). In addition (6), the distance from landmark 1 to landmark 14 (top view) was measured (Fig. 1) to obtain a characteristic variable (SR, standardized radius) instead of the radius, which was not available for all specimens due to the sporadic lack of W4.

Figure 1.

Shell of Pyramidula pusilla s. l. Left: lateral view showing position of landmarks (1–7) and distances measured to gather height of first four whorls (H1–H4). Middle: dorsal view showing position of landmarks (1–14) and distances measured to gather width of protoconch (W0) and first four whorls (W1–W4). Right: ventral view showing landmarks (1–2) used to measure umbilicus diameter (UD).

We used the software tpsUtil v. 1.47 (Rohlf, 2012) and tpsDig v. 2.16 (Rohlf, 2010) to set 13–20 landmarks (cross view, each quarter whorl) from the top view and 5–7 landmarks from the lateral view along the sutures. We used full generalized Procrustes analyses (GPA, data prepared with Kendall coordinates) to quantify shape variation, whereby the shapes are superimposed and size is standardized. All statistical and morphometric analyses were calculated with R v. 2.15.1 (R Development Core Team, 2008). To compare only shape variation, we used baseline registration, in which all landmark configurations (of all specimens) are scaled according to two fixed landmarks, serving as anchors for the superposition, which still enables all possible shapes to be expressed. Thin plate spline (TPS) analyses display the direction and degree of the deviation. For these analyses, we made use of the descriptions and recommendations given by Claude (2008). Multivariate statistics were used to test the affiliation of a specimen to a genetic group and to assess whether there are measurements characteristic for a specific group (logistic regression).

RESULTS

Phylogeographic analysis and mtDNA sequence diversity

The dataset of 357 COI Pyramidula sequences plus one outgroup sequence of Orcula dolium comprised 77 haplotypes. Phylogenetic tree reconstruction (by BI and ML), including the outgroup, resulted in trees with three distinct clades, as well as additional isolated lineages splitting from the basal nodes (data not shown). Because of the low node support, especially for the basal nodes, it was suspected that the outgroup selection was problematic. Orcula dolium had been chosen because of the sister-group relationship between the families Orculidae and Pyramidulidae in the tree presented by Wade, Mordan & Naggs (2006). The assumption that the bad resolution was due to the high distance between outgroup and ingroup sequences (average p-distance 20%) was supported by the fact that the BI and ML trees calculated without outgroup yielded higher support for four distinct clades (1–4). With the exception of clade 3, the relationships between these clades were well supported. Clades 1 and 2 as well as clades 3 and 4 were resolved as sister groups (Fig. 2).

Figure 2.

BI tree based on 76 COI haplotypes of Pyramidula sp. Each lineage is identified by abbreviation of the mountain ranges and number of the specimen. Midpoint rooting was applied. Support values: posterior probability values of the BI analysis (left), bootstrap values of the ML analysis (right).

Dividing the phylogenetic tree into four clades is to some extent arbitrary, since there are two distant haplotypes in clade 1 (Pysp3) and clade 3 (Pypus214) that could be also defined as separate clades. This, however, would have resulted in two clades represented by single individuals. Therefore, we designated only four clades (treating those lineages as part of clades 1 and 3, respectively) to describe the topology and features of the phylogenetic reconstructions as clearly as possible. It is clearly desirable to investigate larger samples from the respective geographic regions of Pysp3 (Greece, Epirus) and Pypus214 (Turkey, Pontic Mountains).

The genetic variation is quite high in all clades, whether measured by p-distances, haplotype diversity (Hd) or nucleotide diversity (π) (Tables 1 and 2). Hd is high in all four groups (0.81–0.86), whereas nucleotide diversity is higher in clades 1 and 3 than in clades 2 and 4. The highest mean distance was observed within clade 1 (2.9%, max. distance 7.0%), whereas the highest mean distances between groups was found between clades 1 and 3 (8.4%, max. distance 9.3%). However, these values are most likely affected by the single distant haplotypes in clades 1 and 3.

Table 1.

Maximum, mean and minimum p-distances (in %) within and between individuals of the four genetic clades of Pyramidula pusilla s. l. and the outgroup (OG) Orcula dolium.

Clade		1			2			3			4		OG
	Max.	Mean	Min.	Max.	Mean	Min.	Max.	Mean	Min.	Max.	Mean	Min.	Mean
1	7.4	2.9	0.0	9.5	7.8	6.7	10.0	8.4	7.1	9.9	7.8	6.5	20.2
2				1.3	0.5	0.0	8.5	7.0	5.2	8.0	7.3	5.7	19.7
3							5.3	1.7	0.0	7.2	5.8	4.2	19.2
4										3.5	0.8	0.0	19.6

Table 2.

Number (n) of sequences and haplotypes, haplotype diversity (Hd) and nucleotide diversity (π) of the four genetic clades and the whole dataset (all) of Pyramidula pusilla s. l.

Clade	1	2	3	4	All
n sequences	14	12	110	221	357
n haplotypes	7	6	21	42	76
Hd	0.81	0.86	0.83	0.84	0.92
π	0.03	0.01	0.02	0.01	0.04

Geographical distribution of the samples

All 14 specimens belonging to clade 1 were collected in the Dinaric Alps (Albania, Greece) and in the Pontic Mountains (Turkey), while clade 2 comprises 12 specimens from the Pontic Mountains and from the Taurus Mountains (Turkey). Clade 3 contains sequences of 110 specimens collected in Austria, Croatia, Germany, Italy, Slovakia, Slovenia, Switzerland and Turkey. Clade 4 includes sequences of 221 specimens from Austria, Germany, Italy, Slovenia and Switzerland (Fig. 3A). Thus, the phylogenetic tree can be roughly divided into two groups: specimens belonging to clades 1 and 2 are without exception from the Balkans or western Asia, whereas specimens belonging to clades 3 and 4 were almost exclusively collected within the Alpine region and surrounding areas, as well as the western Carpathians. Exceptions are two individuals collected in Turkey, one (Pypus202) popping up in the middle of clade 3, the other (Pypus214) splitting from the basal node of clade 3. Repeated DNA extraction, PCR and sequencing of these two samples confirmed their authenticity. Clade 3 includes all specimens from the western Carpathians and many from southern regions, whereas clade 4 comprises mostly specimens from the northeastern and western Alps. However, in general the distributional ranges of clades 3 and 4, including the eastern Alpine samples, overlap to a great extent. The analysis of the phylogenetic tree together with the sampling data revealed that representatives of different clades occur syntopically at 12 sampling sites. This is shown in more detail by the MJ networks of clades 3 and 4, in which the geographical regions are indicated by different colours (Fig. 3B, C).

Figure 3.

A. Sampling sites of Pyramidula. Different coloured dots indicate different genetic affiliation of the specimens collected at each site. Enlargement shows details of sampling sites in Austrian Alps and nearby regions. B, C. Haplotype network of Pyramidula clade 3 (B) and clade 4 (C). The different colours indicate the different geographic regions according to the Alpine Club classification of the eastern Alps (AVE) for Austrian samples and abbreviations of mountain ranges for samples collected outside Austria (see Supplementary Material Table S1 for locality details). Bars on stems indicate single mutations between haplotypes. Size of the circles is proportional to number of individuals sharing same haplotype; small black dots symbolize hypothetic haplotypes not present in dataset. Numbers within/beside the circles indicate the haplotype number.

Morphometric analyses

The question that arises from the results of the genetic analyses is whether specimens clustering in different mitochondrial clades are differentiated morphologically as well. Of the 143 individuals measured, 4 are from clade 1, 5 from clade 2, 44 from clade 3 and 90 from clade 4. The height and width of the fourth whorl could be measured only for a subset of specimens (68 and 53 individuals, respectively). We measured an overall mean SH of 1.76 mm (1.21–2.29 mm, σ = 0.20) with an UD mean of 0.72 mm (0.51–1.02 mm, σ = 0.10) and an average shell radius of 1.08 mm (0.91–1.33 mm, σ = 0.08).

Although some of the measured parameters show distinct arithmetic means in some of the four genetic clades, there is also a high variance of these parameters (Tables 3 and 4). For example, with regard to UD, the specimens belonging to clade 2 show by far the narrowest umbilicus, but the values are still within the ranges found for clades 3 and 4. In addition, they also exhibit the highest shells of all clades, but the ranges of clades 2 and 3 overlap.

Table 3.

Average, minimum and maximum of shell height (SH), umbilicus diameter (UD) and shell radius (SR) of the four genetic clades of Pyramidula pusilla s. l.

	SH	UD	SR
Clade 1	1.62 (1.54–1.66)	0.72 (0.67–0.77)	1.14 (1.06–1.24)
Clade 2	1.92 (1.85–1.99)	0.60 (0.56–0.64)	1.06 (0.89–1.17)
Clade 3	1.65 (1.29–2.18)	0.74 (0.57–1.02)	1.02 (0.78–1.25)
Clade 4	1.81 (1.21–2.29)	0.72 (0.51–1.01)	1.11 (0.80–1.38)

Table 4.

Number of individuals (n), mean (µ), median (x̄), standard deviation (σ), minimum and maximum values of the four whorl heights (H1–H4) and five whorl widths (W0–W4) of the four genetic clades of Pyramidula pusilla s. l.

	Clade	n	µ	$\overline{x}$	σ	Min.	Max.
H1	1	4	0.25	0.26	0.02	0.23	0.27
	2	5	0.30	0.29	0.01	0.29	0.32
	3	44	0.19	0.19	0.03	0.08	0.24
	4	90	0.19	0.18	0.04	0.08	0.26
H2	1	4	0.28	0.27	0.03	0.26	0.32
	2	5	0.34	0.33	0.03	0.32	0.39
	3	44	0.22	0.23	0.04	0.14	0.30
	4	90	0.23	0.24	0.04	0.13	0.31
H3	1	0
	2	0
	3	21	0.44	0.46	0.06	0.32	0.53
	4	47	0.48	0.49	0.06	0.31	0.57
H4	1	0
	2	0
	3	10	0.89	0.89	0.06	0.76	0.99
	4	43	0.94	0.93	0.05	0.86	1.07
W0	1	4	0.23	0.22	0.01	0.22	0.24
	2	5	0.23	0.22	0.03	0.20	0.27
	3	44	0.22	0.22	0.03	0.13	0.27
	4	90	0.22	0.22	0.03	0.13	0.29
W1	1	4	0.23	0.23	0.03	0.21	0.27
	2	5	0.24	0.24	0.02	0.19	0.25
	3	44	0.21	0.20	0.02	0.17	0.26
	4	90	0.23	0.23	0.02	0.16	0.27
W2	1	4	0.28	0.28	0.01	0.27	0.3
	2	5	0.24	0.25	0.03	0.18	0.26
	3	44	0.25	0.25	0.02	0.19	0.31
	4	90	0.28	0.28	0.02	0.21	0.34
W3	1	4	0.40	0.40	0.03	0.37	0.44
	2	5	0.35	0.34	0.03	0.32	0.39
	3	44	0.34	0.34	0.03	0.29	0.40
	4	90	0.38	0.38	0.04	0.30	0.48
W4	1	0
	2	0
	3	18	0.32	0.32	0.04	0.26	0.42
	4	35	0.32	0.32	0.05	0.24	0.40

A principal component analysis (PCA) was performed on the dataset containing all measured variables. The coloured dots in Figure 4A represent individuals from different clades. Even though there are regions of overlap, clade 2 is well separated, clades 3 and 4 are partly separated, while clade 1 is broadly overlapping (mainly) with clade 4. The first two components explain 45% of the total variance. On the same set of individuals a linear discriminant function analysis (LDA) comparing the means of variables was performed to test the differentiation of the four clades. Since it is only possible to include variables that were measured in all individuals, we had to exclude the width and height of the third and height of the fourth whorl from the LDA. Over 91% of the variance can be explained by the first two LDAs. A test of group assignment of the specimens was conducted according to their morphological parameters in the course of the LDA. Based on the variation of those parameters measured in all individuals (SH, UD, H1–H2, W0–W3) 84.6% (121 out of 143) of the specimens were assigned correctly to their respective clade (Fig. 4A). The binomial logistic regression shows that the variables SH, shell radius and UD have the strongest impact on the affiliation. To summarize, even though the clades are not completely separated, according to either the PCA or the LDA, both display a clear (clade 2) or at least partial (clades 3 and 4) separation, whereas clade 1 is not well differentiated.

Figure 4.

A. Left: PCA of shells of Pyramidula showing individuals as coloured dots and core zones (centres of gravity of multidimensional scatter plots) of each clade marked as ellipses. Right: LDA showing each individual as a white symbol with coloured edge; the filled coloured symbols indicate individuals wrongly assigned by the LDA. B. Top left: baseline registration of shells of Pyramidula according to landmarks 5 and 6 (see Fig. 1). Transformity grids with mean shape of each clade (interrupted line) compared with the overall mean shape indicated by coloured dots.

Landmark analyses

The mean centroid sizes (CS) of the lateral landmarks of all four clades were calculated and compared. Clade 1 has the largest CS (814.75) and clade 3 the smallest (713.14). Between these two extremes, clade 2 has a larger CS (807.11) than clade 4 (766.53). The TPS method was employed to display shape variation between representatives of the four clades. We compared the mean shape found in each clade with the overall mean shape and chose to display this variation in form of a deformity grid (Fig. 4B). Again, clade 2 is the most aberrant, whereas clades 3 and 4 show no differentiation at all. A similar picture is found with the baselines analysis (Fig. 4B). The mean shape used for the GPA was estimated from the whole dataset to quantify shape variation by removing all information on size, position and orientation. The GPA showed that clade 2 deviates most from the mean shape calculated from the whole dataset (Supplementary Material Fig. S1).

To assess which parameters are most characteristic of each clade, a logistic regression was calculated (Supplementary Material Fig. S2). According to this, the affiliation to clade 3 declines with increasing SH, SR and W of the shell and H of all whorls, whereas increasing UD increases the probability of belonging to clade 3. The results show that none of the various parameters allows straightforward assignment to any of the clades.

The superimposed data was used to perform a (binomial) logistic regression, with which we tested the assignment of the specimens to a clade based on a combination of measurements (SH, UD, W0–W4) and the altitude (asl). We only tested the assignment to clade 3 and clade 4 and excluded clades 1 and 2, since this statistical model requires similar numbers of specimens. We also performed a multiple logistic regression to test if and how the different variables affect the assignment of the specimens to a group. Since this model gains complexity with the number of variables, we only used the variables that proved to be significant for assignment—UD, SH and SR. We also included altitude, which proved to have no effect on assignment. The probability of belonging to clade 3 rises with increasing SH and SR and with decreasing UD, and vice versa with clade 4.

Altitudinal zonation

Summarizing the genetic and morphometric results, the four clades are differentiated in shape and/or size to some extent. To test whether this finding can be attributed to different altitudinal distributions of clades, possible correlations between measurements and altitude were tested with linear regressions. For this analysis altitudinal zonation after the concept of Kilian, Müller & Starlinger (1993) was applied to the sampling sites. However, this concept is only applicable for samples from Austria or surrounding areas, as it refers especially to regional climate and vegetation. Therefore, 90 individuals were excluded from the analysis, including all specimens from clade 1 and clade 2. The remaining individuals represent the majority (270, 75%) of specimens collected for this study. Sampling sites were assigned to one of the following altitudinal zones: submontane, lower montane, middle montane, high montane, lower subalpine, high subalpine and alpine. The last of these was defined by the authors, because the concept did not apply for sampling sites higher than 2,000 m asl.

However, the linear regression showed that altitude has no effect on the different variables (data not shown). Even though most of the individuals collected at higher altitudes belong to clade 4, there is no significant correlation between the group affiliation and altitude.

DISCUSSION

This study was initially conducted to assess the genetic diversity and intraspecific differentiation of Pyramidula pusilla in the eastern Alps and to elucidate the phylogeographic history of this poorly studied species, especially with respect to Pleistocene refugia and postglacial expansions in this region. The inclusion of specimens from distant areas was intended to provide an estimate of the magnitude of intraspecific divergences. Although clades 1 and 2 (from Turkey, Albania and Greece) are represented by only a few individuals, they provide interesting insights. However, the results of molecular genetic and morphometric analyses show a complex picture and the questions arise whether there might be a second species of Pyramidula occurring in Austria and whether some of the specimens analysed from other countries might also represent distinct species. Despite these suggestions, which requires taxonomic consideration (see below), the results demonstrate the unclear systematics of this genus in general.

Summarizing the results, the phylogenetic analyses based on mitochondrial COI sequences reveal two clades (clades 3 and 4) in Austria, of which one also contains individuals from Turkey. All other samples from the Balkan region (Albania and Greece) and Turkey are well separated in the tree and form two distinct sister clades (clades 1 and 2). There are also morphological differences between clades; although the sample sizes in clades 1 and 2 are small, both differ from all other clades in size as well as shape, whereas clade 3 differs from clade 4 only in size. It must be emphasized that the results have to be corroborated by larger samples.

Whether these clades might represent different species is discussed below. In the first part of the discussion we consider the geographic distribution of clades and haplotypes, with emphasis on the mainly Alpine clades 3 and 4.

Phylogeography and migration patterns

Besides the rough division into two geographical groups (clades 1 and 2 vs clades 3 and 4), no unambiguous geographical pattern can be found, regarding either the distribution of clades or of individuals within clades. This becomes most evident in the network of clade 3 (Fig. 3B), where the most common haplotype is present in specimens from the southeastern and northeastern Alps and the Carpathians. In the haplotype networks of the other clades the picture is similarly mixed up.

For both of the mostly Alpine clades the high genetic diversity within clades is inconsistent with genetic bottlenecks. Thus the assumption that the clades correspond to Pleistocene refugia seems unlikely as one would then expect genetically depleted populations within the formerly glaciated area, or in the peripheral regions from which the recolonization of the Alps might have started. As the networks (especially of clade 3) allow no clear determination of a central ancestral haplotype, an identification of ancestral distributional areas is not possible based on the networks. The present geographic mixture within the Alpine distributional area of clades 3 and 4 suggests that these clades (possibly representing species) might in fact have persisted (and possibly coexisted in some regions) throughout long times over vast regions, with large population sizes and high migration rates. This would explain the high haplotype diversities and intraclade divergences. In addition, the possibility that some populations might have survived in ice-free nunataks within the Alps (Holdhaus, 1954; Gittenberger, 1990; Haase, Esch & Misof, 2013) cannot be ruled out. However, with the present data it cannot be decided whether the high genetic diversity currently found within the eastern Alps can be attributed to persistence of haplotypes in nunataks, in peripheral refugia or in larger and even more distant areas, from which postglacial recolonization and mixing of haplogroups took place. Given the high diversity, much larger samples would be necessary to answer this question.

Despite the uncertainties concerning glacial distribution ranges, a fast postglacial recolonization of large parts of the eastern Alps, which were mostly covered with ice during the last glacial maximum, seems to be likely. In particular, the wide distribution of certain haplotypes is striking (Fig. 3B). This implies that P. pusilla is in fact capable of fast expansion. In general, terrestrial gastropods (Cowie, 1984; Baur, 1993; Kleewein, 1999; Pfenninger & Posada, 2002) are supposed to have very little active dispersal ability. Yet, it is not uncommon for slow and small animals to be transported along a river or water body, to adhere to the plumage or legs of birds or to the fur of animals and thus to be transferred over large distances (Baur, 1993; Gittenberger, 2012). Especially in the case of the very small and lightweight Pyramidula, which is known for its adhesive powers (Klemm, 1974), exozoochory could be a relatively frequent means of long-distance transport. Moreover, intraintestinal transport by predators is also a possible mechanism, as shown by Wada, Kawakami & Chiba (2012) for the land snail Tornatellides boeningi. Another reasonable explanation is wind transportation (Kirchner, Krätzner & Welter-Schultes, 1997). To sum up, frequent passive transport probably played a major role in the recolonization of the Alps by P. pusilla as well as in its biogeographic history in general.

Species delimitation

All specimens analysed from Austria and surrounding countries had initially been assigned to P. pusilla, and those from the Balkan Peninsula and western Asia were tentatively assigned to this species as well. Do the results suggest that the four clades (or some of them) represent distinct species? The high intraspecific distances of 10.96% between the two most distant haplotypes are high at first sight, especially considering that the COI gene is one of the more conserved genes of the mitochondrial genome (Haase et al., 2003). However, high intraspecific distances are not uncommon in terrestrial gastropods, as studies on various species have shown: e.g. around 13% for COI sequences in Arianta arbustorum (Haase et al., 2003), 18% in Orcula dolium (Harl et al., 2014) and 19% in Trochulus hispidus (Kruckenhauser et al., 2014). In the mitochondrial 16S rRNA gene Thomaz, Guiller & Clarke (1996) found genetic distances of 13% in Cepea nemoralis. Several possible reasons for such high genetic distances within snail species have been proposed (e.g. Murray, 1964; Thomaz et al., 1996; Kruckenhauser et al., 2014): high substitution rates, long-lasting isolation followed by secondary contact and reunion of differentiated populations, and finally large population sizes over long periods of time. Pyramidula pusilla fulfils all the requirements for long-lasting persistence of distinct lineages (in our case clades 3 and 4) and consequently higher intraspecific diversity and genetic distances, i.e. a large, fragmented distributional range with a high number of subpopulations together with a presumably low (active) migration rate (as generally assumed for terrestrial gastropods, e.g. Cowie, 1984; Baur, 1993; Kleewein, 1999; Pfenninger & Posada, 2002). Thus, the high genetic distances within and between clades alone are not sufficient to refute the assumption that our complete dataset represents a single species, P. pusilla. However, the fact that the clades are differentiated morphologically and occur sympatrically and even syntopically (at least 1 + 2, 3 + 4) challenges this assumption.

If our sample could represent more than one species, the question arises whether each of the four clades corresponds to a species (i.e. two species within the eastern Alpine region as well as two species within the Balkans or western Asia). The syntopic occurrence of specimens with clade 3 or 4 haplotypes together with the (slight) morphological differentiation can be taken as a hint of the existence of more than one species. But could the differences in size be explained by altitudinal factors? Altitudinal zonation shows that individuals belonging to clade 4 are more likely to be from high (middle montane to alpine) areas than individuals belonging to clade 3. Some authors have observed a correlation between shell size and altitude or altitude-correlated factors such as strong ultraviolet radiation, lower temperatures and shorter seasons (Baur, 1984; Goodfriend, 1986; Baur & Raboud, 1988; Baur & Baur, 1998; Anderson, Weaver & Guralnick, 2007), although others have observed no such correlations (Engelhard & Slik, 1994; Hausdorf, 2003). However, the tests performed in the present study with clades 3 and 4 indicate that the morphological shell characters are independent of altitude, which supports the existence of two species. Moreover, the distribution ranges of clades 3 and 4 overlap widely and there were 11 localities at which both types were collected.

For clades 1 and 2 the sample sizes are low, but the concordance between the gene tree and morphology (size and shape), together with the sympatric occurrence, again argues in favour of two independent species. In general, the concordance between genetics and morphology makes the ‘single species’ hypothesis less likely than the ‘multiple species’ hypothesis. Analyses of ncDNA (e.g. microsatellites) to quantify the extent of gene flow, as well as breeding experiments between different clades, could be used to test these hypotheses.

Taxonomic considerations

If we considered the ‘multiple species’ scenario as more likely, the question arises which clade represents P. pusilla s. s. and which other nominal species could be involved. According to the currently known distributional ranges of the European species of Pyramidula as compiled in Figure 5, P. pusilla is the only species covering the sampling localities in Austria, Germany, Slovakia as well as Turkey. Thus, a reasonable interpretation would be that P. pusilla in fact comprises two cryptic species (represented by clades 3 and 4). An alternative possibility is that the distribution of P. rupestris extends much more to the north than reported so far and that this species is represented by either clade 3 or 4 (P. pusilla would then be the other one).

Figure 5.

Distribution of the six nominal European species of the genus Pyramidula reconstructed after Gittenberger & Bank (1996) and http://www.animalbase.org. Positions of sampling localities are indicated and coloured according to their respective clades.

Concerning clades 1 and 2, the following species occurring on the Balkan Peninsula and in western Asia have to be considered: P. cephalonica, P. pusilla and P. rupestris, as well as P. chorismenostoma. The last of these can be excluded from these considerations because of its conspicuously scalariform shell. Assuming that P. pusilla and P. rupestris correspond to clades 3 and 4, only P. cephalonica remains as a possible name for clades 1 and/or 2. It has a presumed distribution range covering the localities of both clades in the northern parts of Turkey (according to Gittenberger & Bank, 1996). The morphological differentiation (together with cooccurrence) of clades 1 and 2 argue against the assumption that they both represent P. cephalonica. Thus, the question is whether clade 1 or 2 (or both) corresponds to P. cephalonica. However, all these considerations remain speculative, and still other taxonomic outcomes are possible. Of the four genetic clades, specimens belonging to clade 1 have the smallest mean SH (1.62 mm) and the widest mean SR (1.14 mm). According to the morphological descriptions by Gittenberger & Bank (1996), only P. cephalonica and P. umbilicata have a smaller SH than P. pusilla, whereas only P. umbilicata also has a wider SR than P. pusilla. According to these diagnoses, the most likely solution is that specimens represented by clade 1 are P. umbilicata. However, this species has only been described to occur in Great Britain and the northern Iberian Peninsula, raising doubt about this assignment. Concerning clade 2, specimens belonging to this genetic group have the highest mean SH and the narrowest mean UD among the four clades. If we compare our findings with the descriptions by Gittenberger & Bank (1996), only P. jaenensis has a higher SH along with a narrower UD than P. pusilla. Yet, as with clade 1, the currently known geographic distribution of P. jaenensis (southern Iberian Peninsula) would argue against such an assignment. This exemplifies the general problems underlying the taxonomy of Pyramidula: the geographical origin of specimens might be misleading, as the data imply that distributional ranges and possibly even species numbers might have to be revised.

In view of this complexity, it has to be asked whether the original descriptions of species and examination of their types are adequate to permit unambiguous assignment. In general, the original descriptions of the various species are brief and the information on shells (measurements and shape) reported in the literature are not very helpful for differentiation of P. rupestris and P. pusilla. Without information on distributional ranges, assignment based on these descriptions is not possible. Even the type material and type localities do not add significant information. Pyramidula rupestris and P. pusilla were both described rather vaguely and imprecisely localized from ‘France’ (Draparnaud, 1801; Vallot, 1801). The type material of P. pusilla is lost, as reported by Gittenberger & Bank (1996). Examining the type material of P. rupestris at NHMW, we found it to comprise specimens that are very variable in shell morphology, including morphotypes typical of both P. rupestris (i.e. a higher, more slender shell with a narrower umbilicus) and of P. pusilla (a broader and lower shell with a wider umbilicus) (Fig. 6). Measurements of a set of 20 specimens (including lectotype and paralectotypes designated by Gittenberger & Bank, 1996) indicate a comparatively high and slender shell and a narrow umbilicus, concordant with the morphological descriptions of P. rupestris by Gittenberger & Bank (1996) (data not shown). Comparing these data with the measurements of the present investigation, it might appear plausible that clade 4 represents P. rupestris and clade 3 P. pusilla. However, it should be emphasized that both species have been shown here to be very variable in shell morphology, with high overlap in single characters. Thus, the ‘typical’ shapes shown by Gittenberger & Bank (1996) may represent merely two extreme forms out of the range of variation that is found in both clades 3 and 4.

Figure 6.

Lateral and ventral view of paralectotypes of Helix rupestris Draparnaud, 1801 showing that the type material is morphologically variable, from seemingly typical P. rupestris (two individuals on left: NHMW14800, 1820_XXVI_97A_6) to specimens that appear more like P. pusilla (two individuals on right: 1820_XXVI_97A_10, 1820_XXVI_97B_1).

To conclude, our combined genetic and morphological results suggest that the phylogenetic clades might actually represent more than one species and that, contrary to former assumptions, two species of Pyramidula may occur within the eastern Alps. Tentatively we propose that clade 3 represents P. pusilla and clade 4 P. rupestris, although no definitive taxonomic conclusions are yet possible. To resolve the taxonomy of Pyramidula species, analysis of larger samples covering the whole distributional range of the genus, including the described areas and type localities of the various species, will be necessary.

ACKNOWLEDGEMENTS

This project was funded by the Austrian Science Fund (FWF Project No. 19592-B17). We thank the Freunde des Naturhistorischen Museums Wien for financial support for travel expenses and Barna Páll-Gergely, Alexander and Peter Reischütz for providing essential samples. We are grateful to Edmund Gittenberger, Oihana Razkin, Willi Pinsker, Werner Mayer and Lukas Raneburger for valuable discussions and comments on the manuscript and to Hans Winkler for support concerning statistics. Thanks to Barbara Däubl, Julia Schindelar and Oliver Macek for technical assistance in the laboratory.