Phylogeography of the land snail genus Orcula (Orculidae, Stylommatophora) with emphasis on the Eastern Alpine taxa: speciation, hybridization and morphological variation

Background The Central and Southern European mountain ranges represent important biodiversity hotspots and show high levels of endemism. In the land snail genus Orcula Held, 1837 nine species are distributed in the Alps and a few taxa inhabit the Carpathians, the Dinarids and the Western Black Sea region. In order to elucidate the general patterns of temporal and geographic diversification, mitochondrial and nuclear markers were analyzed in all 13 Orcula species. We particularly aimed to clarify whether the Alpine taxa represent a monophyletic group and if the local species diversity is rather the result of isolation in geographically separated Pleistocene glacial refuges or earlier Tertiary and Quaternary palaeogeographic events. In order to test if patterns of molecular genetic and morphological differentiation were congruent and/or if hybridization had occurred, shell morphometric investigations were performed on the Orcula species endemic to the Alps. Results The phylogenetic trees resulting from the analyses of both the mitochondrial (COI, 12S and 16S) and the nuclear (H4/H3) data sets revealed three main groups, which correspond to the three subgenera Orcula, Illyriobanatica and Hausdorfia. The reconstruction of the historic geographic ranges suggested that the genus originated in the Dinarides during the Middle Miocene and first colonized the Alps during the Late Miocene, giving rise to the most diverse subgenus Orcula. Within the latter subgenus (including all Alpine endemics) almost all species were differentiated by both molecular genetic markers and by shell morphometrics, except O. gularis and O. pseudodolium. Conclusions The present study confirms the importance of the Alps as biodiversity hotspot and origin center of land snail diversity. The species diversity in the subgenus Orcula was likely promoted by Miocene to Pliocene palaeogeographic events and the insular distribution of preferred limestone areas. In some cases, speciation events could be linked to the divergence of populations in glacial refuges during the Pleistocene. Sporadic contact between geographically separated and reproductively not yet isolated populations led to intermixture of haplogroups within species and even hybridization and mitochondrial capture between species. Electronic supplementary material The online version of this article (doi:10.1186/s12862-014-0223-y) contains supplementary material, which is available to authorized users.


Background
European mountain ranges harbor a large number of endemic species and are generally considered as important biodiversity hotspots. Concerning diversity in European terrestrial snails, the IUCN Red List of Threatened Species lists 554 native species from rocky areas, 445 species from shrublands and 367 from forest habitats. Among the first category, rocky areas, the Dinarides represent the most diverse European mountain region with about two hundred native gastropods listed, followed by the Alps and the Carpathians, both with somewhat less than a hundred native species listed [1]. Obvious reasons favoring the diversity in mountain areas are the strong structuring of habitats with a wide range of ecological niches, and the availability of different geological substrates. As most land snails are calciphilous, mountain regions offering limestone bedrock are particularly rich in species and show high rates of endemism, whereas intermediate areas with siliciclastic bedrock constitute migration barriers for many taxa. The isolation in favorable habitats is therefore an important reason for diversification of Alpine land snails [2]. Moreover, the current distribution and diversity patterns of Central European land snails and other biota were affected most strongly by climatic events during the Pleistocene. The shifts in temperature and humidity, and the expansion of glaciers, resulted in the fragmentation of populations of many taxa, complete or local extinction, and the loss of genetic variation due to bottlenecks [3]. As large parts of the Alps were covered by glaciers during the Last Glacial Maximum (LGM; 30-18 kya [4]) and earlier glacial periods, their role as origin center of biodiversity was much discussed in the past decades. Reviews of early molecular genetic studies suggest that the Central European mountain ranges did not provide refuges during glacial maxima, but were settled recently from more southern regions [3,5,6]. However, more recent molecular genetic studies support the presence of northern refuges in the periphery of the Alps and in the Western Carpathians for particular organisms like plants and invertebrates [7][8][9][10][11][12]. Populations of former refuge areas were usually characterized by high genetic diversity and the presence of rare (private) alleles [13].
In the present study we investigate the phylogeny and phylogeography of the calciphilous land snail genus Orcula Held, 1837. Orcula species are high-spired snails, 5 to 10 mm in height, with internal lamellae extending to the aperture margins. The morphologies of these lamellae serve as the primary characters for species identification. Currently, 13 species are known, almost all distributed either in the Alps, the Carpathians or in the Dinarides. Only a single species, Orcula zilchi Urbański, 1960, was recorded from the western Black Sea region of Bulgaria and Western Anatolia. Sphyradium dobrogicum Grossu, 1986 was also classified within the genus Orcula in the Fauna Europaea checklist [14], but without any published reference. However, based on available information, the species was synonymized with Sphyradium doliolum (Bruguière, 1792) [15]. Data on type specimens and taxonomic considerations about all Orcula taxa are summarized in the type catalogue of [16].
So far, the most comprehensive investigation of the genus was performed by Gittenberger [17], who attempted to differentiate the Alpine Orcula taxa and the Dinaric Orcula schmidtii (Küster, 1843) by anatomical and shell morphological traits. Schileyko [18] also investigated the genital anatomy and formulated hypotheses about the relatedness of several species. Páll-Gergely et al. [15] were the first to study the anatomy of Orcula jetschini (Kimakowicz, 1883) and O. zilchi and, based on differences in the morphology of the penial caecum, the shell structure and the morphology of the aperture folds, subdivided the genus into three subgenera: (1) Orcula, (2) Illyriobanatica Páll-Gergely & Deli 2013 and (3) Hausdorfia Páll-Gergely & Irikov 2013. The subgenus Orcula includes all species distributed in the Alps, among them the type species Orcula dolium (Draparnaud, 1801), which has by far the widest distribution including the Alps, the Western Carpathians and surrounding lowlands [19]. In contrast, the Alpine endemics are almost exclusively restricted to rocky limestone habitats of the Northern and the Southern Calcareous Alps. Their distribution was mainly investigated by Zimmermann [20] and Klemm [2]. Orcula gularis (Rossmässler, 1837) shows a disjunct distribution in both the Northern Calcareous Alps (Salzburg, Styria and Upper Austria) and the Southern Calcareous Alps (Carinthia and East Tyrol). Orcula austriaca Zimmermann, 1932 shows a similar distribution but its main area is situated more easterly in the Northern Calcareous Alps of Lower Austria. Orcula pseudodolium Wagner, 1912 inhabits a few mountains in the Northern Calcareous Alps of Upper Austria only. A fourth species, Orcula fuchsi Zimmermann, 1931, is restricted to two mountains (Mt. Gippel and Mt. Göller) in the Northern Calcareous Alps of Lower Austria. The other four Alpine endemics are exclusively found in the Southern Calcareous Alps of Austria, Slovenia and Italy. Of these, Orcula tolminensis Wagner, 1912 stands out as it resembles conchologically a dwarf form of O. gularis, with similar aperture characteristics; it is known from three sites in Southern Carinthia and Slovenia only. Orcula restituta (Westerlund, 1887) is mainly distributed in the Slovenian Alps and Orcula spoliata (Rossmässler, 1837) has an isolated distribution about 200 km west in Trentino-Alto Adige (Italy). The fourth Orcula species of the Southern Calcareous Alps is Orcula conica (Draparnaud, 1801). It is common in the eastern part of the Southern Calcareous Alps, but was found at a single site in the Dinarides around the Plitvice lakes (Republic of Croatia) as well. The subgenus Illyriobanatica comprises O. schmidtii, Orcula wagneri Sturany, 1914 and O. jetschini. O. schmidtii and O. wagneri inhabit high mountain regions of the Dinarides, from the Republic of Serbia to southern Greece, and their distribution ranges do greatly overlap [21]. A delimitation of the two taxa is problematic, because the shell characters are highly variable and do not allow to clearly distinguish the two species (see pictures in [16] and [21]). Therefore, the latter two taxa are referred to as O. wagneri/schmidtii complex in the following. O. jetschini was reported from Romania and is the only Orcula species of the Western Romanian Carpathians. It is distributed in the Banat region, western Transylvania (including Crisana) and northern Oltenia. Its shell shows similarities to those of the Dinarid species O. wagneri and O. schmidtii, but it is a woodland species of lower elevations, occurring mainly among leaf-litter or decaying dead wood. Orcula zilchi represents the monotypic subgenus Hausdorfia. It is known from three localities in South-Eastern Bulgaria and from three sites in Western Anatolia only. Its habitat ranges from leaf litter of alluvial forests in the western Black Sea region to limestone rocks in Western Turkey [15,22].
The present paper addresses the evolutionary history of the genus Orcula in general as well as the phylogeographic patterns of the species endemic to the Alps in particular. We aimed to answer the following questions: Which geographic areas were inhabited by ancestral populations of Orcula? What are the causes for the high species diversity in the Alpine region? Is there a congruency between molecular genetic patterns and morphologically defined groups in the Alpine Orcula species? Are there indications for recent or past hybridizations between any of the species?
We performed the first phylogeographic study of the genus Orcula based on comprehensive mitochondrial (mt) and nuclear (nc) data sets including material of all 13 species. A molecular clock analysis was performed and combined with a phylogeographic range reconstruction to trace the distribution patterns of the mt lineages throughout time.
In order to test whether the differentiation in shell morphology is congruent with the molecular genetic groupings, morphometric landmark analyses were conducted with the group of Orcula species endemic to the Alps and the Alpine-Dinarid O. conica.

Phylogenetic trees
A 655 bp fragment of the mitochondrial (mt) cytochrome c oxidase subunit I (COI) was analyzed in 295 specimens from 151 sites ( Figure 1 and Table 1), including samples of all 13 extant Orcula species (Figure 2). The sequences of the Orcula species endemic to the Alps and the Alpine-Dinarid O. conica constituted three quarters of the samples. The nine species were in the focus of the study and we aimed to infer the degree of intraspecific molecular genetic (mtDNA) and morphological variation across the major parts of their distribution areas.
In the Bayesian Inference (BI) and Maximum Likelihood (ML) trees calculated with the COI data set only, most species clades are well supported (data not shown), but the relationships between the clades are not resolved. The complete COI data of the seven Alpine endemics and O. conica are shown separately as phylogenetic networks.
In order to set up a reliable phylogenetic framework, we additionally analyzed sections of the mt 12S ribosomal RNA (12S), 16S ribosomal RNA (16S) and the nuclear (nc) histone H3 and H4 complex (H4/H3) in a selection of 87 individuals. We selected a representative sample of the major mt COI lineages and geographic distributions (including most type localities). Substitution saturation in the mt COI and the trimmed 12S and 16S alignments was examined with the test of Xia et al. [23], implemented in MEGA v.5.1 [24]. The alignments show only little substitution saturation, with I ss.c values significantly larger than I ss values (P = 0.000): 12S (I ss.c 0.699 > I ss 0.218), 16S ( Iss.c 0.722 > Iss 0.249) and COI (I ss.c 0.718 > I ss 0.271). However, moderate substitution saturation (P = 0.0012) is observed in the 3rd codon positions of the COI with I ss.c (0.686) and I ss (0.581) values differing only marginally between each other.
The BI and ML phylograms calculated with the concatenated alignments (COI, 12S and 16S) show congruent, well resolved topologies with high support values for most of the nodes (Figure 3). The topology remains the same when the 3rd codon position of the COI is excluded from the data set, but most nodes show lower posterior probabilities or likelihood values, respectively (Additional file 1). The nc trees were calculated with the H4/H3 data of the same subset of specimens. The overall topology with a division into three main clades representing the three subgenera Orcula, Illyriobanatica and Hausdorfia is illustrated by the nc tree as well ( Figure 4). However, the nc data set is less variable, and support values for most nodes are lower in the nc tree: The first main clade in the mt and nc trees corresponds to the subgenus Orcula, which comprises all Orcula species showing distributions in the Alps. O. dolium is clearly monophyletic and the sister group to a clade comprising the lineages of the Alpine endemics and the Alpine-Dinarid species O. conica. In the mt trees ( Figure 3, Additional file 1), O. conica branches off from the basal node, whereas in the nc tree O. spoliata and O. restituta, the two being sister species, split basally ( Figure 4) (Figure 4), but not in the mt trees ( Figure 3, Additional file 1).

Molecular clock analysis and reconstruction of geographic range history
Molecular clock analyses and reconstructions of the geographic range histories were performed to analyze the temporal and geographic patterns of divergence of the mt lineages. The molecular clock dated trees were calculated  in BEAST v1.7.5 [25] with the three mt markers (COI, 16S and 12S), using Sphyradium doliolum, Orculella bulgarica (Hesse, 1915) and Orculella aragonica (Westerlund, 1897) as outgroups. Since the inclusion of the sequences of the two Orculella species as additional outgroups affected the patterns in the 12S and 16S alignments, the resulting maximum clade credibility trees differ in their topology from the phylograms: The clade of the subgenus Illyriobanatica branches off from the basal node in the molecular clock trees ( Figure 5, Additional files 2 and 3), whereas the subgenus Hausdorfia takes this position in the mt phylograms ( Figure 3, Additional file 1). However, the node marking the first split within the genus Orcula obtained rather low support in all (also the nc histone) trees, and the relation between the three subgenera can still be considered as unresolved.
Three different approaches were performed to estimate the divergence times: (1) In the first approach (Additional file 2), the root of the tree (node XIV) was calibrated to the age of the fossil Nordsieckula falkneri (Hausdorf, 1995), the presumed most recent common ancestor of Orcula, Sphyradium and Orculella. The mean age of the node marking the split of Illyriobanatica from the subgenera Orcula and Hausdorfia (node I) was estimated to 12.2 mya (14.8 to 9.5 mya, 95% HPD interval). The split between the subgenera Orcula and Hausdorfia (node II) was estimated to 9.9 mya (12.6 to 7.1 mya, 95% HPD), and the date of divergence of O. dolium from the other eight Alpine species (node III) was estimated to 8.5 mya (10.9 to 6.0 mya, 95% HPD). (2) In the second approach (Additional file 3), the divergence date of the outgroup taxa O. bulgarica and O. aragonica (node XVI) was calibrated to the time of the first occurrence of ancestral O. aragonica in the fossil record of the Iberian Peninsula. The resulting estimated mean node ages were 13.4 mya (23.5 to 4.7 mya, 95% HPD) for node I, 11.0 mya (19.4 to 3.9 mya, 95% HPD) for node II, and 9.4 mya (16.6 to 3.4 mya, 95% HPD) for node III, respectively. (3) The third approach ( Figure 5) was a combination of the first two and included the calibration of two nodes (XIV and XVI). The estimated mean node ages were 12.4 mya (14.9 to 9.8 mya, 95% HPD) for node I, 10.2 mya (12.9 to 7.5 mya, 95% HPD) for node II, and 8.7 mya (11.2 to 6.3 mya, 95% HPD) for node III, respectively.
The node ages are largely congruent in approaches 1 and 2, suggesting that the two different calibration points did not produce conflicting results. However, the resulting 95% HPD intervals are extremely large when only the split between the outgroup taxa O. bulgarica and O. aragonica (node XVI; approach 2) is calibrated (Additional file 3). In approach 3 (two node datings) ( Figure 5), ranges are similar as in approach 1 (Additional file 2), indicating that the 95% HPD ranges are mainly influenced by placing a prior on the stem of the tree.
The reconstruction of geographic ranges was performed with Lagrange [26], using the linearized maximum clade credibility tree inferred with approach 3 (two dating points). The distribution area of the genus Orcula was classified into seven geographic mountain areas, and the ancestral lineages were allowed to occupy a maximum of two ranges at the same time: In the constrained model (1), migration was prohibited between very distant areas, because some ancestral ranges were biogeographically extremely unlikely, and different migration probabilities were assigned  between adjacent and not immediately adjacent areas ( Figure 5). In the unconstrained model (2), migration was permitted between all areas and with the same dispersal probabilities (Additional file 4).
The results are largely congruent at the outer branches of the trees, but the ancestral ranges and likelihoods estimated for the basal nodes/branches differ strongly between the two approaches. In the constrained model, the most recent common ancestor of the genus Orcula (node I) was likely distributed in the Dinarids (65%), and the ancestor of the subgenera Orcula and Hausdorfia (node II) was most likely distributed in the Dinarids (28%) or in an area additionally including the Bulgarian Strandzha Mts (23%) ( Figure 5). The ancestor of the subgenus Orcula (node III) was most likely distributed in the Western Carpathians and the Dinarids (37%), the Western Carpathians and the Southern Calcareous Alps (13%), or in the Dinarids (11%). In the unconstrained model, a Dinarid ancestry of the genus is the most likely scenario as well, but with a lower probability (25%); an alternative range additionally includes the Northern Calcareous Alps (12%). Accordingly, the common ancestor of the subgenera Orcula and Hausdorfia (node II) was most likely distributed in either the Northern Calcareous Alps (16%) or in the Dinarids (11%). The unconstrained model also predicts different ranges for the common ancestor of the subgenus Orcula (node III): Northern Calcareous Alps (25%), Western Carpathians and the Southern Calcareous Alps (19%), or Western Carpathians and the Dinarids (15%).
The reconstruction of the geographic range history indicates that migrations to geographically distant mountain ranges represented rare events -the Western Black sea area, the Western Alps and the Western Romanian Carpathians were probably colonized only once. The Dinarides were probably re-colonized only once from the Southern Calcareous Alps during the Late Pleistocene or Holocene, namely by O. conica. However, in the subgenus Orcula migrations between Southern and Northern Calcareous Alps seem to have happened several times. Most complex patterns were found in O. dolium, which probably originated in the Western Carpathians and is now found in most limestone areas of the Alps. The results suggest that the species migrated repeatedly between the Western Carpathians and the Alps, and colonized the Western Alps probably only once.

Genetic distances
In the mt data set, distances measured between the subgenera and species are extremely high ( Table 2 Figure). The colored symbols at the branch tips indicate the geographic origin of each haplotype. At the cladogenesis events (nodes), all alternative ancestral subdivision/inheritance scenarios with likelihoods of 10% or more are indicated, separated by an "or", together with the respective likelihoods in%. When scenarios for cladogenesis events involve two ancestral areas, the symbol for the likely ancestral area/-s is/are provided left to each of the two branches. For nodes representing major splits, node ages and 95% posterior HPD intervals are indicated (see table) three Orcula subgenera towards the outgroup S. doliolum are 9.7% (H4: 2.8; H3: 6.4; spacer: 20.0%).

Mitochondrial diversity in the subgenus Orcula
The Orcula species endemic to the Alps and the Alpine-Dinarid O. conica constituted more than three quarters of the samples analyzed and were in the focus of the present study. In order to display distributional patterns of the mt COI haplotypes, Median joining networks were calculated ( Figure 6). The    Figure 6F) shows two highly diverged sequence clusters, separated by 29 mutational steps. The pattern might be the result of a long evolutionary history in the Southern Calcareous Alps.

Morphological variation in the subgenus Orcula
The morphometric analysis was performed with the landmark data of the Alpine Orcula species (except for O. dolium) to evaluate the amount of morphological differentiation between the species and subspecies, respectively. The high intraspecific variability of shell morphs, even from the same localities, complicates a clear separation of the different species based on data of single specimens only. However, a higher resolution is obtained when mean shapes of several specimens per locality are compared with each other (Figure 7). In the Linear Discriminant Analysis (LDA) including the mean shapes of all eight species, 85.02% of the total variance were explained by the first three discriminants (LD1: 59.37%, LD2: 16.09% and LD3: 9.56%; Figure 7A).  Figure 7A). However, when only the landmark data of the latter taxa (and the closely related O. tolminensis) are compared, the species form well-defined clusters ( Figure 7B). In the respective LDA ( Figure 7B), 88.79% of the total variance is explained by the first three discriminants (LD1: 42.33%, LD2: 28.38% and LD3: 18.08%).

Phylogeny and phylogeography of the genus Orcula
One of the main aims of this study was to clarify the phylogenetic relationships among Orcula species and to test whether the Alpine Orcula species represented a monophyletic group. Based on a comprehensive phylogenetic data set including samples of all 13 Orcula species there is clear evidence that the nine species distributed in the Alps represent a monophyletic group corresponding to the subgenus Orcula as proposed by Páll-Gergely et al. [15]. Similarly, Illyriobanatica and Hausdorfia are each monophyletic in the mt and nc trees (Figures 3, 4). Previous considerations about the relationships of the Orcula species were made in particular by Gittenberger [17] and Schileyko [18], based on shell morphological and anatomical traits. Gittenberger's [17] [18] is not consistent with our results and would result in paraphyletic species complexes. The reconstruction of the geographic range history supports a scenario in which the genus Orcula originated in the Dinarides during the Middle Miocene ( Figure 5). Thus, despite the fact that nine of 13 extant Orcula species are distributed in the Alps, the area was most likely not the center of origin, but was colonized from the Dinarids. The lineages of the three subgenera most likely split during the Middle and Late Miocene -during that time period, Alps, Dinarides and Carpathians were partly separated by a lateral branch of the Mediterranean Sea [27]. The separation of O. dolium from the group including the eight Alpine endemics (and the Alpine-Dinarid O. conica) was dated in the Late Miocene and could be explained with the formation of Lake Pannon, which separated Eastern Alps and Western Carpathians during the Tortonian and reached its maximum extent about 10 mya [28]. The results suggest that O. dolium originated in the Western Carpathians and colonized the Eastern Alps first during the Pliocene. The radiation into numerous mt lineages during the Pleistocene can be explained with divergence in separated glacial refuges as suggested by Harl et al. (2014) [19]. The diversification of the other Orcula species endemic to the Alps (including the Alpine-Dinarid O. conica) probably started shortly after the split from O. dolium, during the Late Miocene or the Lower Pliocene. Accordingly, the diversification of most Alpine Orcula species pre-dated the Pleistocene and their speciation cannot plausibly be explained solely by divergence in separate glacial refuges as proposed by Zimmermann (1932) [20] -only the split of the closely related O. austriaca, O. tolminensis and O. gularis is dated in the Pleistocene ( Figure 5).
The species of the subgenus Illyriobanatica are distributed in the Dinarides and the Western Romanian Carpathians. Several other land snail taxa share similar distribution patterns and inhabit both mountain ranges, for instance the hygromiid Xerocampylaea zelebori (L. Pfeiffer, 1853), the aciculid Platyla wilhelmi (A. J. Wagner, 1910) and the clausiliid genus Herilla Adams & Adams, 1855 [29]. Our data suggest that the Dinarid O. wagneri/schmidtii complex and the Western Romanian O. jetschini split during the Late Miocene. The formation of Lake Pannon could have influenced their separation as well -parts of the Western Romanian Carpathians (including the current distribution area of O. jetschini) formed islands during the Middle Tortonian [28,30].
Summarizing, the reconstruction of the geographic range history indicates that the separation of the major groups within the genus Orcula was linked to palaeogeographic events which shaped Europe during the Miocene. The patchy distribution of limestone rock certainly constituted an important factor in the diversification of lineages because all Orcula species are more or less calciphilous. Since lowland areas like the Pannonian Basin featured almost no limestone rock, active migration between mountain ranges was probably hampered. The Pleistocene glaciations obviously had a strong impact on the current distributions of the Alpine endemics because their areas are all located near the eastern margins of the LGM glacier line. However, most of the latter species presumably diverged from each other in the Late Miocene and the Pliocene already, and not during the Pleistocene ( Figure 5).

Hybridization within the subgenus Orcula
In most of the species of the subgenus Orcula, coherent mt and nc sequence patterns (Figures 3, 4, 6) as well as common morphological traits (Figure 7) were found. Thus, these species seem to be reproductively isolated from each other. Moreover, several species were found to occur sympatrically without indication for hybrid-  (Figures 3, 6), and only three single specimens of O. gularis (ENN21_377, ENN23_924 and ENN4_885) from the Ennstal Alps (Austria, Styria) show distinct mt variants, which form a sister clade of O. tolminensis. Specimens from four sites (OOV17, OOV18, OOV19, OOV20) even show transitional states in the expression of the palatal folds (from a backward orientated tooth to a diagonal bulge) (Figure 8), indicating recent hybridization.
Potential causes of non-monophyly of species in mt trees theoretically can be inferred from the depth of the coalescences in gene trees, geographical distribution of shared genetic markers, and concordance with results of admixture analyses of nuclear multilocus markers. [31]. In land snail species, incomplete lineage is discussed in species of the hygromiid genus Xerocrassa Monterosato, 1892 [32] and the helicid Cornu aspersum [33], whereas mt introgression is assumed to have happened in other Xerocrassa species [31] and in the camaenid genus Euhadra Pilsbry, 1890 [34].
The  Figure 3). Shedding more light on this topic would require additional sampling the potential hybridization area and analyzing (additional) nc markers from a larger number of specimens.
Another issue addressed in our study is the morphological variability within O. austriaca. Apart from the common form, three subspecies were described for O. austriaca. Among those, Orcula a. pseudofuchsi is of special interest because its shells are more elongated than those of the nominate form of O. austriaca. Klemm [35] hypothesized that O. a. pseudofuchsi represents an 'intermediate' between O. a. austriaca and O. fuchsi, or descended from the same common ancestor at least. Despite unequivocal shell morphological differences between specimens of O. a. austriaca and O. a. pseudofuchsi (Figure 7), the two taxa could not be delimitated by the mt and nc markers analyzed (Figures 3 and 4). Moreover, no intermediates between O. austriaca and O. fuchsi were found at Mt. Göller (Lower Austria) where both species co-occur and, in contrast to the assumption of Klemm [35], the two species were not even sister species (in the mt trees). Hence, the aberrant shell shape of O. a. pseudofuchsi is most likely not the result of hybridization but rather evolved uniquely in the population of O. austriaca from Mt. Gösing (Lower Austria). Similarly as in O. austriaca, conspecific populations strongly differing in shell morphology but not in their mtDNA were observed in the Western Carpathian populations of O. dolium. Orcula dolium brancsikii Clessin, 1887, exhibiting strongly elongated shells, was found next to populations with specimens featuring rather globulous shells, but a distinction of the two forms was not possible with the nc and mt markers used [19]. Another subspecies, O. a. faueri, inhabits the Southern Calcareous Alps, geographically separated from the nominate form of the Northern Calcareous Alps. Despite the geographic distance, the mt haplotypes of O. a. faueri are embedded within the diversity of the Northern Calcareous Alpine population, suggesting that O. austriaca colonized the Southern Calcareous Alps very recently ( Figure 6). In contrast O. a. faueri could not be differentiated from the common form by its shell morphology (Figure 7). Only in the nc H4/H3 trees there is a comparably strong bifurcation between the lineages of the Northern Calcareous Alps and the Southern Calcareous Alps, indicating that the evolutionary history of the species is probably more complex (Figure 4). A possible explanation could be that the populations of the Northern Calcareous Alps and the Southern Calcareous Alps diverged in allopatry, but intermixed recently, leading to mt capture and the loss of the genuine mt variants of O. a. faueri.

Glacial refuges of the Orcula species endemic to the Alps
Although the diversification of the Orcula species endemic to the Alps probably started already in the Late Miocene and can only partly be attributed to speciation in glacial refuges, the Pleistocene LGM (30-18 kya; [4]) and earlier glacial maxima obviously affected the current distribution to a great extent. Since all of the latter species are strictly calciphilous, they probably could not survive the LGM in lowland areas surrounding the Alps like the closely related O. dolium [19]. However, most of the Alpine endemics show wide altitudinal ranges from the valleys up to high mountain areas (e.g., O. gularis from 400 to 2000 m asl) and are adapted to cold climates. Even though the LGM climatic snowline was 1000 to 1500 meters below the current level in peripheral ranges of the Eastern Alps [36], lower mountain ranges potentially provided suitable conditions during glacial periods. Moreover, the current distributions of all Alpine endemic Orcula species include areas not covered by ice during the LGM (Figure 9). The existence of Eastern Alpine refuges is also supported by several recent molecular genetic studies dealing with mountain plants [8,37,38] and invertebrates [9,10,39]. Besides, the Eastern Alpine margins harbor several endemic species with low active dispersal capabilities, among those blind troglobiotic beetles [40] and endemic land snail species restricted to high altitudes [41]. The genetic diversity within Orcula allows to conclude that the Alpine endemics outlasted the LGM and earlier Pleistocene cold stages probably in several smaller refuges at the periphery of both the Northern and the Southern Calcareous Alps and did not suffer from genetic bottlenecks.

Conclusions
The results of the present study indicate that the evolutionary history of the genus Orcula dates back to the Middle Miocene. The three subgenera most likely derived from an ancestor that was distributed in the Dinarides, and their separation can be explained by palaeogeographic events preventing migration between the mountain ranges populated. Although the Alps were probably not the origin center of the genus, they gave birth to the majority of species. The structuring of the Alps, with two geographically separated major limestone areas (Southern and Northern Calcareous Alps), was of great importance for the diversification of the local Orcula species, most of them being strictly calciphilous. Within the group of Alpine endemics, most speciation events seem to predate the Pleistocene. Their current distribution patterns, however, were strongly shaped by the LGM glacier extent. Most taxa could be differentiated well by both morphologic and genetic traits, with the exception of O. gularis and O. pseudodolium. The latter two species differ in their shell morphology and nc DNA but cannot be distinguished by their mt DNA sequences, indicating mt introgression.

Study area and sampling
Samples of all 13 Orcula species were collected in the years 2007 to 2012 in ten different countries. Figure 1 shows the location of collection sites investigated. The distribution areas in Figure 1 are based on literature data [2,15,21,42], collection data of the Natural History Museum Vienna (NHMW) and the Senckenberg Museum (SMF), and data of the present study. Of a total of 115 localities investigated in the present study, most sites (101) were located in the Alps (Table 1). Elevation and position of the localities were determined via GPS. Samples were collected in various habitats covering an altitudinal range from 190 m to 2200 m above sea level (asl). One to four specimens of each species per site were prepared for the DNA sequence analyses and stored in 80% ethanol, following the protocol of Kruckenhauser et al. [43]. In addition, a number of empty shells were collected from the same sites for the morphometric analyses. DNA samples of the outgroup species Orculella bulgarica and Orculella aragonica were obtained from B. Gómez-Moliner (Universidad del País Vasco, Vitoria, Spain) and were already used by Arrébola et al. [44]. All other voucher specimens were deposited in the Natural History Museum Vienna (NHMW). In order to provide an overview of the taxa investigated, pictures of shells of selected type specimens are shown in Figure 2. The species and subspecies determination was based on a combination of shell characters (expression of aperture folds and shell form) and the geographic distributions reported in literature [2,21].

Outgroup selection for phylogenetic trees and fossil calibration
In the course of the investigations on Orcula, samples of most other orculid genera (including Alvariella Hausdorf, 1996, Orculella Steenberg, 1925, Pilorcula Germain, 1912, Pagodulina Clessin, 1872, Sphyradium Charpentier, 1837 and Schileykula Gittenberger, 1983) were analyzed for the same set of mt and nc markers. Preliminary analyses based on this data set clearly support the monophyly of a group containing Orcula, Schileykula, Sphyradium and Orculella, with Orcula being the sister group to the other three genera (Harl et al. in prep.). Hence, all of these three genera represent equivalent outgroup taxa. We used the monotypic Sphyradium doliolum as outgroup for the calculation of the mt (COI, 12S, 16S) and nc (H4/H3) phylograms (Figures 3, 4). Within the family Orculidae, S. doliolum is by far the most widespread species with an area extending from Western Europe to Kyrgyzstan [45]. Orculella bulgarica and Orculella aragonica were additionally included as outgroups in the molecular clock analyses. O. bulgarica is the most widespread Figure 9 Distribution of mt clades of the Orcula species endemic to the Alps. As the distribution areas of five of the species overlap in the Southern Calcareous Alps, the data are displayed in two separate maps (A and B). The colored symbols indicate the distribution of the mt clades, the color shaded areas represent the distribution areas of the respective Orcula species. In case that two species or species clades, respectively, were found at the same localities, the symbols are shown in ovals pointing towards the locality. The white shaded areas represent the glacier extent during the LGM (30-18 kya; [4]). The data on the maximum extent of glaciers during the LGM was published by [72], and modified by [8].
species within Orculella, distributed from southern Europe to Armenia [45], whereas O. aragonica is the only Orculella species of the Iberian Peninsula [46]. Their sister group relationship is supported by molecular genetic (COI, 16S) and anatomical data [44], and a (within the genus uniquely) shared preference for wet habitats such as small marshes [46]. The paleontological record of O. aragonica comprises more than 30 Pliocene to Holocene sites, of which the Almenara-Casablanca karst complex (Castellón, Spain) features the oldest records (1.8 mya [47]). The site features a continuous Miocene to the Early Pleistocene fossil record, allowing to determine the species' first occurrence in that area fairly precisely. One of the assumptions in the molecular clock analyses was that this first record coincides with the colonization of the Iberian Peninsula by ancestral O. aragonica or the split between O. aragonica and O. bulgarica, respectively. The second assumption was that the fossil orculid species Nordsieckula falkneri represents the most recent common ancestor of Orcula, Orculella, Sphyradium and Schileykula Gittenberger, 1983 (not included in the present study). The species is known from Middle Miocene sediments in Gründlkofen bei Landshut (Germany) and Nowa Wieś Królewska (Poland) and was previously classified into the genus Orcula [48]. Based on strong morphological similarities N. falkneri is considered as closely related to the genera Orcula, Orculella, Sphyradium and Schileykula [1], but as a distinctive trait it features a toothlike subangularis [49], which was yet solely found in some of the extant North-African Orculella species. Its only congeneric Nordsieckula subconica (Sandberger, 1858) is known from Late Oligocene (Hochheim, Germany) to Early Miocene (Tuchořice, Czech Republic) successions [50]. N. subconica features an additional trait, a third columellar lamella (infracolumellaris), which is neither present in N. falkneri nor in any extant species of the genera Orcula, Schileykula, Sphyradium and Orculella ( Figure 10). Since N. falkneri temporally succeeds N. subconica and was found in the same area, it might be a direct descendant of the latter. In the molecular clock analyses, the basal node of the tree was dated with the latest record of N. falkneri (Nowa Wieś Królewska, Poland), which was classified to the Astaracian (MN6/7) [51,52], corresponding to a time period comprising the Middle Miocene Langhian (15.97 to 13.65 mya) and the Serravallian (13.65 to 11.61 mya) [53]. The node age was set to 13.79 mya (SD ±1.0), the mean age between the lower boundary of the Langhian and the upper boundary of the Serravallian.  (4). Unique labels, consisting of a specimen number and a locality tag, were assigned to every specimen. The latter is a three letter code, with each code representing a geographic mountain region (Table 1). DNA was extracted using the QIAgen Blood and Tissue Kit. A section of the mt COI (655 bp) was analyzed in all specimens with the primers COIfolmerFw [12] (modified from Folmer et al. [54]) and H2198-Alb [55]. As the genetic distances between the Orcula species were extremely high, deeper splits in the phylogenetic trees could not be resolved with the COI alone. Therefore, additional mt and nc markers were analyzed in a selection of 86 individuals of Orcula and the eleven outgroup specimens: sections of the mt 12S (673 to 725 bp) and 16S (838-890 bp), and a sequence region comprising large parts of the nc histone genes H3 (345 bp) and H4 (259 bp), and the intermediate non-coding spacer (244-490 bp). The 12S primers (12SGastFw and 12SGastRv) and the H4/H3 reverseprimer (OrcH3Right1) were published in [56]. The 16S primers (16SLOrc1Fw, 16SLOrc2Fw and 16SLOrcRv) and the other H4/H3 primers (OrcH4Left1, OrcH4Left2, OrcH3Sright3 and OrcH4SLeft3) were developed by [19]. For direct sequencing, DNA fragments were amplified with the RocheTaq © DNA polymerase with 3 mM MgCl 2 . The PCR started with 3 min at 94°C, followed by 35 cycles with 30 s at 94°C, 30 s at the particular annealing temperatures (Table 3), 1 min at 72°C, and a final extension for 7 min at 72°C. Primer sequences and annealing temperatures are listed in Table 3. The H4/H3 fragments were amplified with the primers OrcH4Left1 (or the alternative primer OrcH4Left2) and OrcH3Right1, but sequencing was performed with the internal primers OrcH3Sright3 and OrcH4SLeft3. Some specimens showed to be heterozygous regarding the non-coding spacer region. In those cases, the PCR was repeated with the proofreading Finnzymes Phusion © polymerase. The PCR started with 30 s at 98°C, followed by 35 cycles with 10 s at 98°C, 10 s at the particular annealing temperatures (Table 3), 30 s at 72°C, and a final extension for 7 min at 72°C. The PCR products were excised from 1% agarose gels and purified using the QIAquick © Gel Extraction Kit (QIAGEN), extended by A-endings with the DyNAzyme II © DNA polymerase (Finnzymes) and then cloned with the TOPO-TA © cloning kit (Invitrogen). Purification and sequencing (in both directions) was performed at LGC Genomics (Berlin, Germany) using the PCR primers, except for H4/H3 (see above). The primer sequences are shown in Table 2. All sequences are deposited in GenBank under the accession numbers KM188500 -KM188950.

Sequence statistics and phylogenetic tree reconstruction
The raw sequences were edited manually using Bioedit v.7.1.3 [57]. The alignment of the COI sequences was straightforward since there were no insertions or deletions (indels). Median-Joining networks were calculated for the eight Alpine endemics with Network v.4.6.0.0 (Fluxus Technology Ltd.) applying the default settings. In order to reduce unnecessary median vectors the networks were then post-processed with the MP (Maximum parsimony) option.
BI and ML phylograms were calculated based on the concatenated alignments (COI, 12S and 16S) of 86 Orcula specimens and the outgroup S. doliolum. The subset included all single specimens of O. wagneri/schmidtii, O. jetschini, O. zilchi and O. dolium, as well as a selection of specimens from the Alpine endemics. The 12S and 16S sequences were aligned with ClustalX v.2 [58] and adjusted manually. Less conserved sequence regions were excluded by trimming the alignment with TrimAl v.1.3 [59]. The original 12S alignment contained 790 bp of which all 271 gap sites were removed using the 'no gap' option (removal of all sites containing gaps). Another 67 positions were excluded by applying the 'strict' option (trimming based on an automatically selected sequence similarity threshold [59]). The original 16S alignment contained 945 positions of which 241 gap sites were excluded. Another 111 sites were removed with the 'strict' option. Subsequently, COI, 12S and 16S were concatenated and identical sequences were collapsed, resulting in a total of 86 unique haplotypes.
Substitution saturation was assessed for all single mt data sets using the test of [23], implemented in DAMBE v.5.2.78 [60]. In order to accommodate substitution saturation in the 3rd codon positions of the COI and to test the influence on the phylogeny, alternative trees were calculated excluding the third codon positions of the COI. Having collapsed identical sequences, the concatenated alignment of this data set contained a total of 77 haplotypes.
The optimal substitution models were determined for all individual data sets with JModelTest v.2.1.5 [61], based on the corrected Akaike Information Criterion (AICc). TrN + I + G was the best-fit substitution model for COI, COI '1st and 2nd codon positions' and 16S, and TPM1uf + I + G for 12S. The optimal model for the entire concatenated alignments was GTR + G + I. However, owing to the limited number of models applicable in MrBayes [62,63], the evolutionary model was set to GTR + G + I for all separate data partitions in the Bayesian analyses.
BI and ML phylograms were also calculated with the H4/H3 data set. Similarly, the alignment was split into three partitions, namely H4, spacer and H3. The spacer was aligned with ClustalX v.2 [58] and all sites with gaps were removed. The data set contained 94 sequences (including additional clones when specimens were heterozygous for the H4/H3) of 87 specimens, which were collapsed to 53 unique haplotypes. Based on the AICc, the best fitting substitution models were K80 for H3 (345 bp) and H4 (259 bp), and K80 + G for the noncoding spacer (206 bp, gaps excluded). The optimal model for the concatenated H4/H3 sequences was GTR + G.
The BI analyses were calculated using the concatenated (mt and nc) alignments, with three data partitions each, allowing MrBayes v.3.2.2 [62,63] to evaluate the model priors of each partition independently. Applying the respective model parameters, the analyses were run for 5x10 6 generations each (2 runs each with 4 chains, one of which was heated), sampling every hundredth tree. Tracer v.1.5 [25] was used to assess whether the two runs had converged and when the stationary phase was reached, which was the case already after several thousand generations. In a conservative approach, the first 25% of trees were discarded as burnin and a 50% majority rule consensus tree was calculated from the remaining 37,500 trees.
ML bootstrap trees were calculated with MEGA v.5.1 [24], applying the models GTR + G + I to the mt and GTR + G to the nc sequences, respectively, but without using separate data partitions (since this option is not supported by MEGA v.5.1). For all data sets, 500 bootstrap replicates were performed using Subtree-Pruning-Regrafting (SPR) as heuristic method for tree inference.
Based on the alignment including all COI sequences, calculations of mean p-distances between species clades and maximum p-distances within the clades were performed with MEGA v.5.1 [24]. For the 12S and 16S alignments (gap sites excluded), mean p-distances between the species clades were calculated, and for the H4/H3 alignments mean p-distances were calculated between the subgenera only. Haplotype and nucleotide diversities were evaluated with DnaSP v.5.10 [64] for the complete COI data set.
The sequence alignments used for the calculation of the mt and nc trees are provided in Additional files 8 and 9.

Molecular clock analysis and reconstruction of geographic range history
The calculations of molecular clock dated linearized BI trees were performed in BEAST v.1.7.5 [25] with the concatenated mt sequences of COI, 12S and 16S. Apart from S. doliolum (1 specimen), 4 specimens of O. bulgarica and 6 specimens of O. aragonica were included as outgroup taxa, because these were used for dating the trees. The molecular clock analyses were calculated with the complete COI and the trimmed 12S and 16S alignments, following the same procedure as for the inference of the phylograms. Having removed all gap sites and performing the "strict" option in TrimAl v.1.3 [59], the 12S and 16S alignments contained 449 and 575 positions, respectively. For model selection and molecular clock analyses, identical sequences were collapsed to a total of 85 haplotypes (out of 97 specimens). The best fitting substitution models were calculated separately for each partition with JModeltest 2.1.5 [61], based on the AICc, resulting in the models HKY + G + I for COI and TN93 + G + I for12S and 16S. The relative rate variation among lineages was tested separately for all three partitions with the molecular clock test implemented in MEGA v.5.1 [24] applying the optimal substitution models inferred with JModeltest 2.1.5. Since the null hypothesis of equal evolutionary rates throughout the trees were rejected at a 5% significance level (P =0), divergence times were estimated under a relaxed molecular-clock in all molecular clock analyses. The divergence times of the mitochondrial lineages were estimated using three different approaches: (1) The basal node of the tree was dated to 13.79 ± 1.0 (SD) mya (15.43 to 12.15 mya, 95% HPD interval), corresponding to the mean age between the lower boundary of the Langhian (15.97 mya) and the upper boundary of the Serravallian (11.61 mya) [53] or the Middle Miocene Astaracien, respectively, the period to which the fossils of N. falkneri were dated (see chapter 'Outgroup selection for phylogenetic trees and fossil calibration') (Additional file 2). (2) The node marking the split between O. bulgarica and O. aragonica was dated with the presumed first appearance of ancestral O. aragonica in the Iberian Peninsula around 1.8 mya as reported by [47]. The node age was set to 1.8 ± 0.1 (SD) mya (1.96 to 1.64 mya, 95% HPD interval) (Additional file 3). (3) A combination of approaches 1 and 2, applying the dating of nodes mentioned above ( Figure 5).
In the sites settings of BEAUti v.1.7.5 (part of the BEAST package [25]), the best fitting substitution models were applied separately to each of the three partitions and the node datings were assigned in the prior settings. The speciation model Yule Process [65] was chosen as tree prior. The BEAST analyses were each performed with four independent runs for 10 7 generations and every thousandth tree was sampled. Having checked whether the four runs had converged, these were combined with LogCombiner v.1.7.5 (part of the BEAST package). Subsequently, 25% of the trees were discarded as burnin and maximum clade credibility trees were calculated each from the remaining 30,000 trees.
The reconstruction of the geographic range history was performed in Lagrange [26], which uses a dispersalextinction-cladogenesis (DEC) modeling for analyzing ML probabilities of rate transitions as a function of time. The rate-calibrated linearized tree from the molecular clock analysis (approach 3, two dating points) was prepared for Lagrange configurator [26], together with a range matrix in which each taxon/lineage was assigned to one of seven geographic regions, Northern Calcareous Alps (N), Southern Calcareous Alps (S), Western Alps (W), Western Carpathians (C), Dinarides (D), Western Romanian Carpathians (R) and Strandzha Mts. (Z). The maximum number of areas allowed for ancestors (= areas inhabited at the same time) was set to '2'. We tested two different models: (1) High migration probabilities ('1.0') were assigned to migration between directly adjacent areas, and lower probabilities ('0.5') were assigned to migration between not immediately adjacent areas (W to C; N to D; N to R). Migration was prohibited between geographically very distant areas (W to D; W to R; Z to W, N, S and C) ( Figure 5). (2) Migration was permitted between all seven geographic regions with the same dispersal probabilities ('1.0') (Additional file 4).
The sequence alignments used for the calculation of the molecular clock tree are provided in Additional file 10.

Morphological analyses
A total of 526 specimens were analyzed in the morphometric analyses. The sample included only specimens of the group of Alpine endemics: O. austriaca