- Research article
- Open Access
Phylogeography of the land snail genus Orcula(Orculidae, Stylommatophora) with emphasis on the Eastern Alpine taxa: speciation, hybridization and morphological variation
BMC Evolutionary Biologyvolume 14, Article number: 223 (2014)
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.
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.
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.
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 . 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 . 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 . As large parts of the Alps were covered by glaciers during the Last Glacial Maximum (LGM; 30-18 kya ) 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 ,,. 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 -. Populations of former refuge areas were usually characterized by high genetic diversity and the presence of rare (private) alleles .
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 , but without any published reference. However, based on available information, the species was synonymized with Sphyradium doliolum (Bruguière,, 1792) . Data on type specimens and taxonomic considerations about all Orcula taxa are summarized in the type catalogue of .
So far, the most comprehensive investigation of the genus was performed by Gittenberger , who attempted to differentiate the Alpine Orcula taxa and the Dinaric Orcula schmidtii (Küster, 1843) by anatomical and shell morphological traits. Schileyko  also investigated the genital anatomy and formulated hypotheses about the relatedness of several species. Páll-Gergely et al.  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 . 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  and Klemm . 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 . 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  and ). 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 ,.
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.
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. , implemented in MEGA v.5.1 . The alignments show only little substitution saturation, with Iss.c values significantly larger than Iss values (P = 0.000): 12S (Iss.c 0.699 > Iss 0.218), 16S (Iss.c 0.722 > Iss 0.249) and COI (Iss.c 0.718 > Iss 0.271). However, moderate substitution saturation (P = 0.0012) is observed in the 3rd codon positions of the COI with Iss.c (0.686) and Iss (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). In the mt trees, the most basal nodes are weakly supported, while the relationships between the other Alpine endemics are well resolved: O. fuchsi, endemic to the Northern Calcareous Alps, is the sister group of a highly supported clade comprising O. austriaca, O. gularis, O. pseudodolium and O. tolminensis. Within this clade, O. gularis is paraphyletic because most specimens of O. gularis and O. pseudodolium show up in the same clade. Moreover, three specimens of O. gularis (IndIDs 377, 885 and 924) from three distinct sites (ENN4, ENN21 and ENN23; Styria, Austria) possess different mt haplotypes, being closely related to those of O. tolminensis. O. austriaca forms the sister clade of the latter two clades. Differing in a few substitutions or indels only, the H4/H3 sequences do not allow resolving the relationships within the group of the latter species (Figure 4). However, this part of the H4/H3 tree shows some peculiar patterns. There is a rather deep split between samples of O. austriaca from the Northern and the Southern Calcareous Alps, whereas these geographically isolated populations are barely differentiated in the mt sequences. There is also a larger diversity of nc haplotypes within O. pseudodolium when compared to the more widespread O. gularis.
The second main clade corresponds to the subgenus Illyriobanatica and includes O. schmidtii and O. wagneri from the Dinarides, and O. jetschini from the Western Romanian Carpathians. The sequences of O. schmidtii and O. wagneri form a highly supported clade in both the mt and the nc trees, but the two species are not monophyletic. Consequently, they are referred to as O. schmidtii/wagneri complex in the present study. O. jetschini is clearly separated from the O. schmidtii/wagneri clade, but the monophyly of the subgenus Illyriobanatica is well supported (Figures 3, 4). Apart from other patterns of sequence similarity, all taxa of the subgenus Illyriobanatica lack a section of approximately 230 bp in the non-coding spacer region of the H4/H3 sequences.
A third major lineage is constituted by O. zilchi, representing the monotypic subgenus Hausdorfia. It inhabits the western Black Sea region and is not closely related to any other species in the trees. Similarly as the taxa of the subgenus Illyriobanatica, O. zilchi features a very short branch in the H4/H3 tree (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  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 , 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.
In the mt data set, distances measured between the subgenera and species are extremely high (Table 2; Additional files 5 and 6). The uncorrected mean p-distances between the subgenera Orcula and Illyriobanatica are 28.6 (12S), 22.0 (16S) and 21.1% (COI), respectively. The mean distances between Orcula and Hausdorfia are 28.1 (12S), 25.2 (16S) and 24.7% (COI), whereas mean distances between Illyriobanatica and Hausdorfia are 31.4 (12S), 25.6 (16S) and 24.0% (COI), respectively. Surprisingly, the latter distances are even higher than the average distances between the genus Orcula and the outgroup S. doliolum, which are 30.9 (12S), 25.6 (16S) and 22.1% (COI). The species providing the largest intraspecific distances in the mt genes is O. dolium with 15.6 (12S), 14.0 (16S) and 18.3% (COI), followed by the O. wagneri/schmidtii complex with 12.3 (12S), 11.6 (16S) and 15.5% (COI). Within the group of Alpine endemics, highest intraspecific distances are observed in O. tolminensis with 3.7 (12S), 4.8 (16S) and 4.9% (COI), in the clade comprising O. gularis and O. pseudodolium with 4.3 (12S), 3.9 (16S) and 4.6% (COI), and in O. austriaca with 1.5 (12S), 1.9 (16S) and 3.6% (COI), respectively. Haplotype and nucleotide diversities calculated for the separate species clades with the complete COI data set are high in all groups (Table 3). The sequence divergences are also high within the set of H4/H3 sequences (Additional file 7). The mean p-distances are 4.3% (H4: 1.6; H3: 2.1; spacer: 8.5%) between Orcula and Illyriobanatica, and 3.7% (H4: 1.4; H3: 2.4; spacer: 7.2%) between Orcula and Hausdorfia. Contrary to the pattern in the combined mt trees, the branch lengths of Illyriobanatica and Hausdorfia are shorter in the tree calculated with the nc sequences (Figure 4), resulting in a considerably lower mean distance of 2.9% (H4: 1.6; H3: 2.1; spacer: 5.1%) between the two subgenera. Mean distances of the 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 COI network corresponding to the O. gularis/O. pseudodolium clade in the mt trees (Figure 3, Additional file 1) comprises the majority of sequences (132) (Figure 6A1). Three specimens of O. gularis from the Northern Calcareous Alps (Ennstal Alps), clustering with O. tolminensis in the mt trees, are shown as a separate small network (Figure 6A2) because they are too diverged. Alternative scenarios explaining the relation of O. gularis and O. pseudodolium are discussed in detail in the following section. The main network of O. gularis/O. pseudodolium (Figure 6A1) is roughly divided into two clusters of haplotypes: the first includes samples of O. gularis from the Southern Calcareous Alps (Karawanks) and the Northern Calcareous Alps (Ennstal Alps, Dachstein, Ybbstal Alps and Totes Gebirge), and of a single specimen of O. pseudodolium (Northern Calcareous Alps: Upper Austrian Prealps). The second cluster includes all samples of O. pseudodolium (Northern Calcareous Alps: Upper Austrian Prealps), as well as a few haplotypes of O. gularis specimens from both the Northern (Salzkammergut Mts., Tennen Mts. and Totes Gebirge) and the Southern Calcareous Alps (Karawanks). The existence of Southern Calcareous Alpine haplotypes in both clusters indicates at least two independent migration events between the Southern Calcareous Alps and the Northern Calcareous Alpine populations. Since the Southern Calcareous Alpine haplotypes in the second cluster are highly derived, it can be assumed that O. gularis was already present in that area for a longer period, probably before the LGM (30-18 kya; ). However, the similarity of the Southern and Northern Calcareous Alpine haplotypes in the first cluster implies a second, recent migration event from the latter area to the Southern Calcareous Alps. A second clade, shown in a separate network, is formed by three specimens of O. gularis from the Northern Calcareous Alps (Ennstal Alps).
A further network (Figure 6B) comprises all 56 samples of O. austriaca, including sequences of the three subspecies O. a. pseudofuchsi Klemm, 1967, O. a. faueri Klemm, 1967 and O. a. goelleri Gittenberger, 1978. The haplotypes of O. a. faueri (Southern Calcareous Alps: Karawanks) are embedded within a larger variety of haplotypes of the Northern Calcareous Alpine populations of O. austriaca. They are separated only by a single mutational step from haplotypes present in the Gutenstein Alps. The haplotypes of O. a pseudofuchsi (Northern Calcareous Alps: Gutenstein Alps) are located in two different clusters of the network, one with haplotypes of the Southern Calcareous Alpine O. a. faueri, the second with haplotypes of a neighboring Northern Calcareous Alpine mountain area (Rax-Schneeberg). The haplotypes of O. a. goelleri (Mürzsteg Alps) also show up in two different clusters of the network. Concluding, none of the three subspecies of O. austriaca is clearly delimited from the nominate form in its mt COI sequences.
The networks representing O. fuchsi (Figure 6E) and the Southern Calcareous Alpine endemics O. conica (Figure 6D), O. spoliata and O. restituta (Figure 6C) are less complex. The sequences of the latter two species are shown in a combined network because we found only one haplotype in O. spoliata. Noticeable is that the Dinarid (Mala Kapela, Croatia) specimens of O. conica feature a haplotype which was found also in the Southern Calcareous Alps (Karawanks), more than 200 km south-east, indicating a recent long distance migration event. The network of O. tolminensis (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). In the LDA, shell shapes of O. conica, O. fuchsi and O. tolminensis are clearly differentiated. A distinct cluster is formed by O. restituta and O. spoliata, whose shells strongly resemble each other. The shell shapes of O. a. pseudofuchsi are unique and differ clearly from those of the nominate form O. a. austriaca, while O. a. faueri is hardly differentiated. The shape clusters of O. gularis, O. austriaca and O. pseudodolium are slightly overlapping (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. . 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  and Schileyko , based on shell morphological and anatomical traits. Gittenberger’s  suggestion that O. gularis, O. tolminensis, O. pseudodolium and O. austriaca are close relatives, demarcated from species representing rather independent lineages (O. dolium, O. spoliata, O. restituta, O. conica and O. fuchsi), was confirmed in the present study. In contrast, the phylogenetic scheme proposed by Schileyko  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 . 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 . 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) . 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)  - 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 . 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 ,.
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 hybridizations: O. austriaca and O. tolminensis, O. fuchsi and O. austriaca, O. gularis and O. conica, O. dolium with O. austriaca, and O. gularis with O. conica and O. fuchsi. Nonetheless, our data strongly support that hybridization happened between O. gularis and O. pseudodolium. Although these two species could be discriminated by their shell forms in the morphometric analyses (Figures 7A and 7B) and provided different H4/H3 sequences (except for a single specimen of O. pseudodolium from the potential hybridization area) (Figure 4), a clear assignment to one or the other species was not possible based on the mt sequences. Most mt haplotypes of O. gularis cluster with O. pseudodolium (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. . In land snail species, incomplete lineage is discussed in species of the hygromiid genus Xerocrassa Monterosato, 1892  and the helicid Cornu aspersum, whereas mt introgression is assumed to have happened in other Xerocrassa species  and in the camaenid genus Euhadra Pilsbry, 1890 .
The non-monophyly of O. gularis and O. pseudodolium in the mt trees might be explained by three different scenarios at least: (1) O. gularis acquired its mitochondria from O. pseudodolium by mt introgression, but genuine mt variants of O. gularis still exist in the population of the Ennstal Alps. (2) O. pseudodolium acquired its mitochondria from O. gularis by mt introgression, and the aberrant mt variants in the three O. gularis specimens of the Ennstal Alps were acquired from O. tolminensis by mt introgression. (3) The mixed mt patterns in the mt O. gularis/O. pseudodolium clade resulted from incomplete lineage sorting. Scenario 1 would require hybridization between O. gularis and O. pseudodolium only and, thus, provide a more parsimonious explanation than scenario 2. A close relationship between O. austriaca, O. gularis and O. tolminensis is also supported by similarities in the genital anatomy, whereas O. pseudodolium shows quite distinct traits . Furthermore, O. gularis and O. tolminensis strongly resemble each other in their aperture traits by the unique presence of a palatal tooth. Incomplete lineage sorting (scenario 3) would not explain the relation between O. pseudodolium and O. gularis sufficiently. The mt sequence patterns in the network (Figure 6A) rather indicate that hybridization occurred at different points of time - some haplotypes of O. gularis are highly diverged from those of O. pseudodolium whereas others are identical or differ only by a few substitutions from each other. Moreover, the existence of the second O. gularis clade in the Ennstal Alps virtually cannot be explained by incomplete lineage sorting (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  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 , 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 . 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 Orculaspecies 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; ) 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. 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 , 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 ,, and invertebrates ,,. Besides, the Eastern Alpine margins harbor several endemic species with low active dispersal capabilities, among those blind troglobiotic beetles  and endemic land snail species restricted to high altitudes . 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.
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 ,,,, 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. . 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. . 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 ,.
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 . Orculella bulgarica and Orculella aragonica were additionally included as outgroups in the molecular clock analyses. O. bulgarica is the most widespread species within Orculella, distributed from southern Europe to Armenia , whereas O. aragonica is the only Orculella species of the Iberian Peninsula . Their sister group relationship is supported by molecular genetic (COI, 16S) and anatomical data , and a (within the genus uniquely) shared preference for wet habitats such as small marshes . 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 ). 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. Based on strong morphological similarities N. falkneri is considered as closely related to the genera Orcula, Orculella, Sphyradium and Schileykula, but as a distinctive trait it features a tooth-like subangularis , 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 . 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) ,, 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) . 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.
PCR and sequencing
A total of 295 specimens were analyzed by means of molecular genetics (specimen numbers in brackets): O. austriaca (56), O. conica (19), O. dolium (35), O. fuchsi (11), O. gularis (84), O. jetschini (2), O. pseudodolium (48), O. restituta (8), O. spoliata (4), O. tolminensis (6), O. wagneri/schmidtii (9), O. zilchi (3), S. doliolum (1), O. bulgarica (6) and O. aragonica (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  (modified from Folmer et al. ) and H2198-Alb . 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 reverse-primer (OrcH3Right1) were published in . The 16S primers (16SLOrc1Fw, 16SLOrc2Fw and 16SLOrcRv) and the other H4/H3 primers (OrcH4Left1, OrcH4Left2, OrcH3Sright3 and OrcH4SLeft3) were developed by . For direct sequencing, DNA fragments were amplified with the RocheTaq© DNA polymerase with 3 mm MgCl2. 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 . 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.220.127.116.11 (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  and adjusted manually. Less conserved sequence regions were excluded by trimming the alignment with TrimAl v.1.3 . 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 ). 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 , implemented in DAMBE v.5.2.78 . 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 , 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 ,, 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  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 non-coding 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 , to evaluate the model priors of each partition independently. Applying the respective model parameters, the analyses were run for 5x106 generations each (2 runs each with 4 chains, one of which was heated), sampling every hundredth tree. Tracer v.1.5  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 , 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 . 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  for the complete COI data set.
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  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 , 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 , 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  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)  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 . 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 ), 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  was chosen as tree prior. The BEAST analyses were each performed with four independent runs for 107 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 , which uses a dispersal-extinction-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 , 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.
A total of 526 specimens were analyzed in the morphometric analyses. The sample included only specimens of the group of Alpine endemics: O. austriaca (66 specimens/12 localities), O. a. faueri Klemm, 1967 (60/4), O. a. pseudofuchsi Klemm, 1967 (21/2), O. conica (22/4), O. fuchsi (16/3), O. gularis (143/19), O. pseudodolium (141/14), O. restituta (18/2), O. spoliata (5/1) and O. tolminensis (34/4). Photographs of the shells were taken in frontal position with a WILD MAKROSKOP M420 and a NIKON DS Camera Control Unit DS-L2, all with the same magnification. Using tpsUtil v.1.44 , tps-files were created, then landmarks were set on 17 uniquely defined points of the shells’ projections in tpsDig v.2.12 . Landmarks 1, 2 and 9 (Figure 7D) are type I, landmarks 16 and 17 type III, and all other points are type II landmarks as defined in . Qualitative characters such as the expression of the aperture folds are not considered in the morphometric analyses. Since the shell forms differ quite among specimens of the same sample localities and, thus, increase the variability of the data set, consensus shapes (here: mean shapes of several specimens of a single locality) were created with tpsSuper v.1.14 . Linear discriminant analyses (LDA) were performed with the landmark data of all Alpine endemics. Plots of the mean shapes were calculated with R v.3.2.0  by using RStudio v.0.97.551 . Morphometric analyses of O. dolium are not performed here, this topic will be discussed in detail in a forthcoming publication (Harl et al. in prep).
Availability of supporting data
The data sets supporting the results of this article are included within the article and its supplementary files.
IUCN Red List of Threatened Species Version 2013.2.., [http://www.iucnredlist.org/]
Klemm W: Die Verbreitung der rezenten Land-Gehäuse-Schnecken in Österreich.. Denkschr Sot Akad Wiss Math-Nat Kl, Volume 117. 1974, 1-503.
Hewitt G: Genetic consequences of climatic oscillations in the Quaternary. Phil Trans R Soc London Ser B Biol Sci. 2004, 359: 183-195. 10.1098/rstb.2003.1388.
Ivy-Ochs S, Kerschner H, Reuther A, Preusser F, Heine K, Maisch M, Kubik PW, Schlüchter C: Chronology of the last glacial cycle in the European Alps. J Quat Sci. 2008, 23: 559-573. 10.1002/jqs.1202.
Hewitt G: The genetic legacy of the Quaternary ice ages. Nature. 2000, 405: 907-913. 10.1038/35016000.
Taberlet P, Fumagalli L, Wust-Saucy A-G, Cosson J-F: Comparative phylogeography and postglacial colonization routes in Europe. Mol Ecol. 1998, 7: 453-464. 10.1046/j.1365-294x.1998.00289.x.
Benke M, Brändle M, Albrecht C, Wilke T: Pleistocene phylogeography and phylogenetic concordance in cold-adapted spring snails (Bythinella spp.). Mol Ecol. 2009, 18: 890-903. 10.1111/j.1365-294X.2008.04073.x.
Schönswetter P, Stehlik I, Holderegger P, Tribsch A: Molecular evidence for glacial refugia of mountain plants in the European Alps. Mol Ecol. 2005, 14: 3547-3555. 10.1111/j.1365-294X.2005.02683.x.
Dépraz A, Hausser J, Pfenninger M, Cordellier M: Postglacial recolonization at a snail’s pace (Trochulus villosus): confronting competing refugia hypotheses using model selection. Mol Ecol. 2008, 17: 2449-2462. 10.1111/j.1365-294X.2008.03760.x.
Weigand AM, Pfenninger M, Jochum A, Klussmann-Kolb A: Alpine crossroads or origin of genetic diversity? Comparative phylogeography of two sympatric microgastropod species. PLoS One. 2012, 7: e37089-10.1371/journal.pone.0037089.
Scheel BM, Hausdorf B: Survival and differentiation of subspecies of the land snail Charpentieria itala in mountain refuges in the Southern Alps. Mol Ecol. 2012, 21: 3794-3808. 10.1111/j.1365-294X.2012.05649.x.
Duda M, Sattmann H, Haring E, Bartel D, Winkler H, Harl J, Kruckenhauser L: Genetic differentiation and shell morphology of Trochulus oreinos (Wagner, 1915) and T. hispidus (Linnaeus, 1758) (Pulmonata: Hygromiidae) in the northeastern Alps. J Moll Stud. 2011, 77: 30-40. 10.1093/mollus/eyq037.
Provan J, Bennett KD: Phylogeographic insights into cryptic glacial refugia. Trends Ecol Evol. 2008, 23: 564-571. 10.1016/j.tree.2008.06.010.
Fauna Europaea version 2.6.., [http://www.faunaeur.org/]
Pall-Gergely B, Deli T, Irikov A, Harl J: Subgeneric division of the genus Orcula Held 1837 with remarks on Romanian orculid data (Gastropoda, Pulmonata, Orculidae). Zookeys. 2013, 301: 25-49. 10.3897/zookeys.301.5304.
Harl J, Sattmann H, Schileyko A: Types of the extant taxa of the landsnail genus Orcula Held 1837 (Gastropoda: Stylommatophora: Orculidae). Arch Molluskenkd. 2011, 140: 175-199.
Gittenberger E: Beitr ge zur Kenntnis der Pupillacea. VIII. Einiges über Orculidae. Zool Verh. 1978, 163: 1-44.
Schileyko A: On the anatomy of Orculidae with special reference to the spermatophores (Gastropoda Pulmonata, Stylommatophora). Ruthenica. 2012, 22: 141-158.
Harl J, Duda M, Kruckenhauser L, Sattmann H, Haring E: In search of glacial refuges of the land snail Orcula dolium (Pulmonata, Orculidae) - An integrative approach using DNA sequence and fossil data. PLoS One. 2014, 9: e96012-10.1371/journal.pone.0096012.
Zimmermann S: über die Verbreitung und die Formen des Genus Orcula Held in den Ostalpen. Arch Naturgeschichte. 1932, 1: 1-57.
Hausdorf B: Zum Vorkommen der Gattung Orcula Held in Griechenland (Gastropoda: Orculidae). Arch Molluskenkd. 1987, 118: 51-55.
Páll-Gergely B: New and little-known landsnails from Turkey (Gastropoda: Pulmonata). Zool Middle East. 2010, 50: 89-94. 10.1080/09397140.2010.10638416.
Xia X, Xie Z: DAMBE: Data analysis in molecular biology and evolution. J Heredity. 2001, 92: 371-373. 10.1093/jhered/92.4.371.
Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S: MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol. 2011, 10: 2731-2739. 10.1093/molbev/msr121.
Drummond AJ, Rambaut A: BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007, 7: 214-10.1186/1471-2148-7-214.
Ree RH, Smith SA: Maximum likelihood inference of geographic range evolution by dispersal, local extinction, and cladogenesis. Syst Biol. 2008, 57: 4-14. 10.1080/10635150701883881.
Harzhauser M, Piller WE: Benchmark data of a changing sea - Palaeogeography, Palaeobiogeography and events in the Central Paratethys during the Miocene. Palaeogeogr Palaeoclimatol Palaeoecol. 2007, 253: 8-31. 10.1016/j.palaeo.2007.03.031.
Harzhauser M, Latal C, Piller WE: The stable isotope archive of Lake Pannon as a mirror of Late Miocene climate change. Palaeogeogr Palaeoclimatol Palaeoecol. 2007, 249: 335-350. 10.1016/j.palaeo.2007.02.006.
Deli T, Subai P: Revision der Vitrea-Arten der Südkarpaten Rumäniens mitBeschreibung einer neuen Art (Gastropoda, Pulmonata, Pristilomatidae). Contrib Nat Hist. 2011, 19: 1-53.
Rögl F, Steininger FF: Vom Zerfall der Tethys zu Mediterran und Paratethys. Die neogene Paläogeographie und Palinspastik des zirkum-mediterranen Raumes. In Ann des Naturhistorischen Museums Wien Ser A für Mineral und Petrogr Geol und Paläontologie, Anthropol und Prähistorie, 1981:135–163.
Sauer J, Hausdorf B: Reconstructing the evolutionary history of the radiation of the land snail genus Xerocrassa on Crete based on mitochondrial sequences and AFLP markers. BMC Evol Biol. 2010, 10: 299-10.1186/1471-2148-10-299.
Sauer J, Hausdorf B: Sexual selection is involved in speciation in a land snail radiation on Crete. Evol (N Y). 2009, 63: 2535-2546.
Guiller A, Madec L: Historical biogeography of the land snail Cornu aspersum: a new scenario inferred from haplotype distribution in the Western Mediterranean basin. BMC Evol Biol. 2010, 10: 18-10.1186/1471-2148-10-18.
Shimizu Ueshima R: Y: Historical biogeography and interspecific mtDNA introgression in Euhadra peliomphala (the Japanese land snail). Heredity (Edinb). 2000, 85: 84-96. 10.1046/j.1365-2540.2000.00730.x.
Klemm W: Über Ostalpine Orculae. Arch Molluskenkd. 1967, 96: 101-111.
Tribsch A, Schönswetter P: Patterns of endemism and comparative phylogeography confirm palaeoenvironmental evidence for Pleistocene refugia in the Eastern Alps. Taxon. 2003, 52: 477-497. 10.2307/3647447.
Parisod C, Besnard G: Glacial in situ survival in the Western Alps and polytopic autopolyploidy in Biscutella laevigata L. (Brassicaceae). Mol Ecol. 2007, 16: 2755-2767. 10.1111/j.1365-294X.2007.03315.x.
Schorr G, Pearman PB, Guisan A, Kadereit JW: Combining palaeodistribution modelling and phylogeographical approaches for identifying glacial refugia in Alpine Primula. J Biogeogr. 2013, 40: 1947-1960.
Homburg K, Drees C, Gossner MM, Rakosy L, Vrezec A, Assmann T: Multiple glacial refugia of the low-dispersal ground beetle Carabus irregularis: Molecular data support predictions of species distribution models. PLoS One. 2013, 8: e61185-10.1371/journal.pone.0061185.
Holdhaus K: Die Spuren der Eiszeit in der Tierwelt Europas. Abhandlungen Zool Gesellschaft Wien. 1954, 18: 1-493.
Duda M, Kruckenhauser L, Haring E, Sattmann H: Habitat requirements of the pulmonate land snails Trochulus oreinos oreinos and Cylindrus obtusus endemic to the Northern Calcareous Alps, Austria. Ecomont. 2010, 2: 5-12. 10.1553/eco.mont-2-2s5.
Lisický MJ: Mollusca Slovenska (The Slovak Molluscs). 1991, VEDA vydavateľstvo Slovenskej akad mie vied, Bratislava
Kruckenhauser L, Harl J, Sattmann H: Optimized drowning procedures for pulmonate land snails allowing subsequent DNA analysis and anatomical dissections. Ann Naturhist Mus Wien B. 2011, 112: 173-175.
Arrébola J, Razkin O, Gómez-Moliner B, Páll-Gergely B: Redescription ofOrculella aragonica(Westerlund 1897), an Iberian species different fromO. bulgarica(Hesse 1915) (Gastropoda: Pulmonata: Orculidae).North West J Zool 2012, 8:.,
Hausdorf B: Die Orculidae Asiens (Gastropoda: Stylommatophora). Arch Molluskenkd. 1996, 125: 1-86.
Garrido JA, Arrebola JR, Bertrand M: Extant populations of Orculella bulgarica (Hesse, 1915) in Iberia. J Conchol. 2005, 38: 653-
Robles F, Martínez-Ortí A: Moluscos continentales de los alrededores de Molina de Aragón (Guadalajara, España), con notas sobre Orculella bulgarica (Hesse, 1915) (Gastropoda, Orculidae). Iberus. 2009, 27: 99-105.
Hausdorf B: Eine miozäne Orcula aus Mitteleuropa (Gastropoda: Orculidae). Heldia. 1995, 2: 73-74.
Harzhauser M, Neubauer TA, Georgopoulou E, Harl J: The Early Miocene (Burdigalian) mollusc fauna of the North Bohemian Lake (Most Basin). Bull Geosci. 2014, 89: 819-908. 10.3140/bull.geosci.1503.
Wenz W: Gastropoda extramarina tertiaria. Fossilium Catalogus I: Animalia. 1923, W. Junk, Berlin, 737-1068.
Nordsieck H: Fossile Clausilien, V. Neue Taxa neogener europäischer Clausilien. II Arch Molluskenkd. 1981, 111: 63-95.
Nordsieck H: Zur Stratigraphie der neogenen Fundstellen der Clausiliidae und Triptychiidae Mittel-und Westeuropas (Stylommatophora, Gastropoda). Mitt Bayer Staatsslg Palaont hist Geol. 1982, 22: 137-155.
Gibbard PL, Head MJ, Walker MJC: Formal ratification of the Quaternary System/Period and the Pleistocene Series/Epoch with a base at 2.58 ma. J Quat Sci. 2010, 25: 96-102. 10.1002/jqs.1338.
Folmer O, Black M, Heah W, Lutz R, Vrijenhoek R: DNA primers for amplification of mitochondrial cytochrome C oxidase subunit I from diverse metazoan invertebrates. Mol Mar Biol Biotechnol. 1994, 3: 294-299.
Gittenberger E, Piel WH, Groenenberg DSJ: The Pleistocene glaciations and the evolutionary history of the polytypic snail species Arianta arbustorum (Gastropoda, Pulmonata, Helicidae). Mol Phylogenet Evol. 2004, 30: 64-73. 10.1016/S1055-7903(03)00182-9.
Cadahía L, Harl J, Duda M, Sattmann H, Kruckenhauser L, Fehér Z, Zopp L, Haring E: New data on the phylogeny of Ariantinae (Pulmonata, Helicidae) and the systematic position of Cylindrus obtusus based on nuclear and mitochondrial DNA marker sequences. J Zool Syst Evol Res. 2014, 52: 163-169. 10.1111/jzs.12044.
Hall TA: BioEdit: a user-friendly biological sequences alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symp Ser. 1999, 41: 95-98.
Larkin MA, Blackshields G: ClustalW and ClustalX version 2. Bioinformatics. 2007, 23: 2947-2948. 10.1093/bioinformatics/btm404.
Capella-Gutiérrez S, Silla-Martínez JM, Gabaldón T: trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics. 2009, 25: 1972-1973. 10.1093/bioinformatics/btp348.
Xia X, Xie Z: DAMBE: data analysis in molecular biology and evolution (version 18.104.22.168). J Hered. 2001, 92: 371-373. 10.1093/jhered/92.4.371.
Darriba D, Taboada GL, Doallo R, Posada D: jModelTest 2: more models, new heuristics and parallel computing. Nat Meth. 2012, 9: 772-10.1038/nmeth.2109.
Huelsenbeck JP, Ronquist F: MRBAYES: Bayesian inference for phylogeny. Bioinformatics. 2001, 17: 754-755. 10.1093/bioinformatics/17.8.754.
Ronquist F, Huelsenbeck JP: MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003, 19: 1572-1574. 10.1093/bioinformatics/btg180.
Librado P, Rozas J: DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009, 25: 1451-1452. 10.1093/bioinformatics/btp187.
Gernhard T: The conditioned reconstructed process. J Theor Biol. 2008, 253: 769-778. 10.1016/j.jtbi.2008.04.005.
Rohlf F: tpsUtil. Version 1.52. 2012. ., [http://life.bio.sunysb.edu/ee/rohlf/software.html]
Rohlf F: tpsDig2. Version 2.16. 2010. ., [http://life.bio.sunysb.edu/ee/rohlf/software.html]
Bookstein F: Morphometric Tools for Landmark Data. Geometry and Biology. 1991, Cambridge Univ Press UK, Cambridge
Rohlf F: tpsSuper. Version 1.14. 2004. ., [http://life.bio.sunysb.edu/ee/rohlf/software.html]
R Core Team: R: A Language and environment for statistical computing v.3.2.0. 2013.
R Studio: RStudio: Integrated development environment for R v.0.97.551. In 2012.
Van Husen D: LGM and late-glacial fluctuations in the Eastern Alps. Quat Int. 1997, 38-39: 109-118. 10.1016/S1040-6182(96)00017-1.
This project was financed by the Austrian Science Fund (FWF Proj.-No. 19592-B17) J. Harl was also supported by the Austrian "Theodor Körner Fonds" (2012) and the University of Vienna (PhD Completion Grant 2013). The ‘Freunde des Naturhistorischen Museums Wien’ provided financial support for travel expenses. The study was also supported by scholarships from the Japan Student Services Organization and the Mitsubishi Corporation. We thank Simone Latkolik and Barbara Däubl for technical assistance in the lab. Many thanks to Benjamín Gómez-Moliner, Tamás Deli, Atanas Irikov and Rajko Slapnik who provided essential samples. Alexander Reischütz, Peter Reischütz, Judith Harl, Florian Mittendorfer, Werner Mayer and Wilhelm Pinsker contributed valuable information and critical comments on the manuscript. Finally, we would like to thank the three anonymous reviewers for valuable suggestions and critical comments.
The authors declare that they have no competing interests.
JH conceived the study, carried out the molecular genetic studies and sequence analyses, created the graphics and drafted the manuscript. BP-G provided essential samples and participated in drafting the manuscript. SK performed the statistic part of the morphological analyses. MD participated in the design and coordination of the study. LK participated in the design and coordination of the study. HS conceived the study and participated in its design and coordination. EH conceived the study, participated in its design and coordination, and essentially helped to draft the manuscript. All authors participated in the field work and read and approved the final manuscript.
Electronic supplementary material
Additional file 1: BI tree of the concatenated mitochondrial sequences ( 12S , 16S and COI 1st2nd positions). Posterior probabilities and ML bootstrap values are provided for all nodes above species level. The scale bar indicates the expected number of substitutions per site according to the models of sequence evolution applied. The black dots indicate nodes with high BI posterior probabilities (1.0) and ML bootstrap values (≥95). The colors of the species clades/labels correspond to those used in Figures 4, 7 and 9. (JPEG 965 KB)
Additional file 2: Linearized molecular clock dated tree (approach 1). Maximum-clade-credibility tree calculated with BEAST, using the concatenated alignments of 12S, 16S and COI. The root of the tree (node XIV) was calibrated to the age of the fossil Nordsieckula falkneri, the presumed most recent common ancestor of Orcula, Sphyradium and Orculella. Node bars indicate the 95% HPD ranges estimated for each node. Mean node ages and 95% HPD ranges are provided for the major splits. A time scale in mya is given below the tree. Abbreviations of the geological epochs: C: Chattian, A: Aquitanian, B: Burdigalian, L: Langhian, S: Serravallian, T: Tortonian, M: Messinian, Z: Zanclean, P: Piacenzian, G: Gelasian and C: Calabrian. (JPEG 1004 KB)
Additional file 3: Linearized molecular clock dated tree (approach 1) . Maximum-clade-credibility tree calculated with BEAST, using the concatenated alignments of 12S, 16S and COI. The divergence date of the outgroups O. bulgarica and O. aragonica (node XVI) was calibrated to the time of the first occurrence of ancestral O. aragonica in the fossil record. Node bars indicate the 95% HPD ranges estimated for each node. Mean node ages and 95% HPD ranges are provided for the major splits. A time scale in mya is given below the tree. Abbreviations of the geological epochs: C: Chattian, A: Aquitanian, B: Burdigalian, L: Langhian, S: Serravallian, T: Tortonian, M: Messinian, Z: Zanclean, P: Piacenzian, G: Gelasian and C: Calabrian. (JPEG 1 MB)
Additional file 4: Reconstruction of the historic geographic ranges (unconstrained model). The linearized molecular clock dated maximum-clade-credibility tree shows the relationships of selected mt lineages (concatenated 16S, 12S and COI sequences). Migration was permitted between all areas and with the same dispersal probabilities. Black dots indicate nodes with high posterior probabilities. 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). A time scale in mya is given below the tree. Abbreviations of the geological epochs: C: Chattian, A: Aquitanian, B: Burdigalian, L: Langhian, S: Serravallian, T: Tortonian, M: Messinian, Z: Zanclean, P: Piacenzian, G: Gelasian and C: Calabrian. (JPEG 1 MB)
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.