Cold tolerant species often exhibit specific responses to Pleistocene climate oscillations, including range expansions during glacial periods and altitudinal shifts between glacial and interglacial periods. Small terrestrial mammals are ideal candidates to study these processes due to their limited dispersal abilities, short generation time and rapid sequence evolution. The aim of this study was to investigate phylogeographical structure of the Alpine shrew (Sorex alpinus) within its recently fragmented range in the central and southern European mountain system. To examine its genetic structure, we sequenced the hypervariable domain of the mitochondrial control region in 51 individuals from 18 localities, covering the majority of the distributional area of the species. We analyzed the sequence dataset using population and landscape genetic approaches. We discovered shallow mitochondrial genealogy with degree of differentiation among site-specific haplogroups. This pattern, together with demographic analyses showing population expansion, corroborate the hypothesis that during the glacial periods, S. alpinus expanded its range into lower elevations between northern and high mountain system glaciers. This expansion was followed by a postglacial range breakup, retreat to moderate elevations of high mountain systems, and a formation of relict enclaves in smaller hilly areas.
Introduction
Cold-adapted arctic and alpine species often exhibit inverse Pleistocene dynamics compared to temperate taxa, usually expanding their ranges during glacial periods (Schmitt 2007, Stewart et al. 2010). Many of these species, known as glacial relicts, survived the warming climate at the end of the last glacial period in fragmented populations restricted to interglacial refugia with characteristic arcto-alpine disjunct ranges (Bhagwat & Willis 2008). Analogical pattern was discovered in cold-tolerant temperate taxa, which might survive glacial periods in higher latitudes within small patches of climatically favorable conditions termed northern glacial refugia (Stewart & Lister 2001, Kotlík et al. 2006, Schmitt 2007). Persistence of temperate species of small mammals in northern refugia is facilitated by their biological traits, such as short generation time, high metabolic rate, cold-tolerance, and habitat generalism (Bhagwat & Willis 2008). Due to the dry nature of the glacial climate, many of these species were restricted to the wetter regions near the foothills of glaciated mountain systems during the glacial periods (Schmitt 2007, Hope et al. 2011), which promoted the frequent origin of different genetic lineages in individual European mountain systems. Populations of above-mentioned species might undergo episodes of altitudinal shifts or fragmentation in response to the post-glacial warming climate, rather than long-distance expansion (Schmitt 2007, Bhagwat & Willis 2008). There is biogeographic evidence for several links between particular mountain systems in Europe, which served as a bridge for the dispersion of animal and plant species (Schmitt 2007).
Many lineages in the orders Erinaceomorpha and Soricomorpha were particularly influenced by+ Pleistocene climatic oscillations, chiefly due to their high metabolic rate and subsequent trophic requirements for high food intake, resulting in dependence on ecosystems with high productivity. However, the character of adaptive response to climate change depends upon dispersal capacities of the respective group. For example, hedgehogs from the genus Erinaceus underwent pulses of range contractions and expansions and became a model group for the role of large Mediterranean refuges in the phylogeography of Europe (Seddon et al. 2001, Bolfíková & Hulva 2012). On the other hand, Mediterranean populations of small mammals like Sorex shrews often show substantial degree of endemism and low contribution to genetic diversity of the species outside the Mediterranean refugia (Bilton et al. 1998). Considering present common occurrence of the genus Sorex in temperate and boreal regions, often crossing the polar circle, these shrews are good candidates to endure glacial period also in northern refugia (Vega et al. 2010).
The Alpine shrew (Sorex alpinus Schinz, 1837) is a small soricomorph with disjunct range in the Alps, the NW Balkans, the Carpathians, and some isolated mountain areas in Germany, Czech Republic and Poland (Spitzenberger 1990, Meinig 2004). It has been probably extirpated from the Harz Mts. (Gahsche 1994) and Pyrenees (Mitchell-Jones et al. 1999). It is a terrestrial and semifossorial shrew restricted to altitudes between 200 and 2550 m above sea level (Spitzenberger 1990). It occupies a variety of humid habitats, usually along mountain streams, in humid crevices of rocky areas and stone walls and also in caves (e.g. Kratochvíl & Grulich 1950, Hutterer 1982, Spitzenberger 2001). Its recent distribution is interpreted as relict (Spitzenberger 1990) and presumably results from range dynamics during the Pleistocene (Hewitt 2004).
The Alpine shrew holds the basal position on the phylogenetic tree of the subgenus Sorex (Fumagalli et al. 1999, Ohdachi et al. 2006), which is of Palaearctic origin (Hutterer 2005). Its fossil history is not well known, but the fossil records, originating mainly from Late Pleistocene, come from within its current range in Bosnia and Herzegovina, Croatia, Slovenia, Germany and Italy (Rzebik-Kowalska 1998). The oldest evidence from the Middle Pleistocene of northern Italy suggests the Alpine origin (Spitzenberger 2001). Classical taxonomy recognized three subspecies of the Alpine shrew described according to the morphological characters: the nominotypic one, S. alpinus alpinus Schinz, 1837 described from the Alps (St. Gotthard-Pass, Suisse); S. alpinus hercynicus Miller, 1909 described from specimens from Harz Mts. and Krkonoše Mts. (i.e. Sudetic Mts.) and S. alpimis tatricus Kratochvíl & Rosický, 1952. The authors designated the population from the High Tatra Mts. as terra typica of the later subspecies and stated that it also inhabits several other Slovak and Czech mountains belonging mostly to the Carpathian mountain system (Kratochvíl & Rosický 1952). However, Spitzenberger (1990), reviewed extensive material from the entire distribution of the species and found great overlap in body- and skull measurements between individual populations. Accordingly, she suggested not to consider S. a. tatricus as a separate subspecies.
Shrews are good adepts for revealing effects of climate oscillation in Pleistocene due to their rapid sequence evolution, which is due to their high metabolic rate and short generation time (Stewart & Baker 1994). In this study we aimed to examine the phylogeographic structure of the Alpine shrew, which is endemic to mountains of Central Europe. We hypothesize that the recent disjunct range is a relict of a wider glacial distribution. Since S. alpinus is dependent on cold and humid conditions, we may use its response to glacial-interglacial cycles as a proxy to address the history of mountains and uplands habitats and associated fauna in Europe.
Material and Methods
Sampling
Tissue samples for this study were collected from 2012 to 2014 in the range of S. alpinus. Tissue samples were stored in 96 % ethanol until DNA extraction. The samples comprised 71 individuals from 18 localities, which covers the majority of S. alpinus distribution (Fig. 1). Elaborated material includes Sudetic Mts. (Krkonoše Mts. and Jizerské hory Mts.), the Carpathians, eastern Alps, the Dinaric Mts. (Mt. Snežnik, Šator Mts., Bjelašnica Mts., Zelengora Mts. and Komovi Mts.) and two isolated small ranges: Žd'árské vrchy Mts. situated between Sudetic Mts. and Šumava Mts., and Rhön Mts. in SE Germany (Fig. 1).
Fig. 1.
Extant range of the Alpine shrew (modified from www.iucnredlist.org) depicted by violet with sampling localities shown by dots.
DNA extraction and sequencing
Genomic DNA was extracted from tissue samples using the commercial DNeasy® Blood and Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. For PCR amplification of the left domain of mitochondrial control region (Fumagalli et al. 1996), primers L15926 (5′TCAAAGCTTACACCAGTCTTGTAAACC3′) from Kocher et al. (1989) and H16498 (5′CCTGAAGTAAGAACCAGATG3′) from Southern et al. (1988) were used. The 25 µl PCR reaction mixtures contained 1× Taq buffer, 2.5 mM MgCl2, 200 µM dNTPs, 0.5 µM primers, 1 U Taq polymerase (Promega, Madison, WI, U.S.A.) and approximately 50–100 ng of genomic DNA. The temperature profile for 35 cycles of PCR amplification was 3 min at 95 °C (initial denaturation step), 1 min at 95 °C, 1 min at 55 °C (annealing temperature), 1 min at 72 °C, followed by a final extension step of 4 min at 72 °C. The reaction was purified using the QIAquick® PCR Purification Kit (Qiagen) or excised from agarose gel and purified with the QIAquick® Gel Extraction Kit (Qiagen). The purified products were sequenced with forward and reverse primers on the 3130 Genetic Analyzer by Applied Biosystems® using the BigDye® Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems®, Waltham, MA, U.S.A.).
Table 1.
Details of sampling localities and number of sampled individuals.
Data analysis
Amplification and sequencing reactions of mitochondrial control region were successful in 71 individuals of the Alpine shrew. Sequences were compiled in the software Geneious®, version 6.1.8. ( http://www.geneious.com, Kearse et al. 2012). An alignment of the sequences was performed using the online program MAFFT, version 7 (Katoh & Standley 2013), which automatically selects an appropriate algorithm according to the data size (Katoh et al. 2002). From the dataset of 71 samples, 20 samples were excluded due to the poor-quality signal of the PCR product or due to the presence of nuclear copies of the mtDNA (numts) in the targeted DNA fragment. Numts were distinguished from the mtDNA sequences as they formed a clade distant from the rest of the sequences (p-distance ranged from 0.094 to 0.109). The numts also contained indels and lacked geographic variation among localities.
High degree of sequence conservation in d-loop-like nuclear copies, contrasting with rapid mutation rate in mitochondrial d-loop, was ascertained also in other studies (Zhang et al. 2006). Consequently, 51 sequences were used for further phylogeographic analyses as final dataset. Three regions with low consensus values resulting from a large amount of missing data or high degree of sequence polymorphisms consisting of indels were excluded to avoid incorrect identification of homologies. Two regions with low consensus values (56 bp at the 5′ end and 106 bp at the 3′) were excluded because of a large number of missing sequences. One section of 316 bp from the R1 variable region (Fumagalli et al. 1996) was excluded due to occurrence of indels (e.g. Warnow 2012). Alignment (both unedited and edited versions) is accessible at TreeBASE ( http://purl.org/phylo/treebase/phylows/study/TB2:S18503). Locality descriptions and sample sizes are given in Table 1. Sequences were collapsed to haplotypes in DnaSp version 5.10.1 (Librado & Rozas 2009) and the haplotypes were submitted to GenBank (accession numbers: KT725725–KT725775).
Relationships among mitochondrial haplotypes were visualized using a median-joining network implemented in the software NETWORK, version 4.613 (Bandelt et al. 1999). Nucleotide polymorphism and haplotype diversity was calculated in the program DnaSP, version 5.10.1 (Librado & Rozas 2009). Best model of nucleotide evolution of ingroup was evaluated in JModelTest 2.1.7 (Darriba et al. 2012). Bayesian inference analysis implemented in the program MrBayes (Ronquist et al. 2012) was used to reconstruct phylogenetic relationships among haplotypes. Analysis consisted of two runs with four chains (one cold and three heated) and was run for 20000000 generations of Markov Chain Monte Carlo (MCMC) simulations sampled every 100 generations. Bum-in period was set to 25 %. Convergence of the two runs was verified in AWTY (Wilgenbusch et al. 2004). Runs were checked for stationarity and convergence by the cumulative plots that show posterior probabilities of splits at selected increments over an MCMC run, and by the bivariate plot of the split frequencies. Resulting Bayesian tree was edited in program FigTree v1.4.2. ( http://tree.bio.ed.ac.uk/software/figtree/).
Landscape genetic method implemented in R (R Core Team 2015) package Geneland (Guillot et al. 2008) was used to assess the spatial distribution of genetic variability of haplotypes in the dataset. Initially, the analysis was performed five times with 100000 MCMC iterations, a thinning of 100, number of populations (K) varying from K = 1 to K = 10, and with an uncorrelated allele frequency model. Finally, the analysis was re-run for the best K with 1000000 iterations, a burn-in of 100 iterations and a thinning of 1000, other parameters had identical settings. Geneland recognized four groups in our data and these groups were treated separately for calculations of the summary statistics of genetic variability which were computed in DnaSP version 5.10.1 (Librado & Rozas 2009).
The occurrence of isolation by distance (IBD) was tested by calculating correlation between matrices of p-distances and matrices of geographic distances among localities. Matrix of p-distances between sequences was calculated in MEGA 5 (Tamura et al. 2011) and matrices of geographic distances among individuals were estimated in Geographic distance matrix generator 1.2.3. (Ersts 2011). Regression line was calculated using the RMA method (Reduced Major Axis). The significance of the correlation was tested with a Mantel test on Isolation by distance web service V3.23 (Jensen et al. 2005). IBD was tested by plotting the parameter p-distance against the logarithm of geographic distance.
The Bayesian skyline plot (BSP) which was implemented in the software Beast v1.8.2. (Drummond et al. 2012) was used to reveal past changes in population size for the entire dataset as well as for each group obtained in landscape genetic analysis separately. The HKY + G model was selected in jModelTest 2.1.7 (Darriba et al. 2012). Analysis of BSP was conducted with 30000000 iterations of the MCMC simulations three times, strict molecular clock with mutation rate of 1, genealogy and parameters of the model were stored every 1000 iterations and burn-in was set to 10000000 iterations. The convergence of chains was assessed in the program Tracer v1.6.0. (Rambaut et al. 2014). We also estimated changes in population size in the past using the mismatch distribution test implemented in DnaSP. Mismatch distribution test compares the observed frequencies of pairwise differences among sequences with the frequencies expected under demographic models. Statistical evaluation of mismatch analysis was conducted in program Arlequin 3.5.2 (Excoffier et al. 2005).
To further assess whether the species experience past population growth we used summary statistics of genetic variability of sequences, including neutrality tests (Fu and Li's F*, Fu and Li's D*, Fu's F's and Tajima's D), computed in DnaSP and Arlequin.
Fig. 2.
Bayesian tree based on 41 haplotypes of Sorex alpinus. Colours correspond to the results of the tesselation map from Geneland (Fig. 4). Numbers at nodes indicate posterior probabilities. Sample codes correspond to those in Table 1.
Fig. 3.
The median-joining network of 41 haplotypes of Sorex alpinus. The size of the circle is proportional to the frequency of the haplotype in the dataset, colour codes represent the geographic origin. Small red dots represent hypothetical haplotypes not found in our sample. Haplotypes from isolated populations at Žd'árské vrchy (CZ-KR) and Rhön Mts. (DE-RH1, DE-RH2), and the CZ-OD5 haplotype (shared by four individuals) are labeled. Numbers at branches represent numbers of mutational steps (displayed for n > 1).
Fig. 4.
Tesselation map illustrating spatial distribution of genetic variability for K = 4, inferred by software Geneland. Black dots represent sample sites. The first (white) group comprises populations from the north-western part of the species range (Sudetic Mts.). Second group (green) comprises populations occurring in the southwestern part of the area (Eastern Alps and Northern Dinarids). The third group (yellow) consists of populations from the southern Dinarids. Fourth group (orange) contains populations from the east part of the species range (the Carpathians).
Results
The comparison of 51 mitochondrial sequences with the length 579 bp revealed 41 different haplotypes with 85 polymorphic sites and overall nucleotide diversity of 0.0199. Genetic distance (p-distance) between individuals ranged from 0.5 % to 4.2 %. Best substitution model was HKY + G according to Akaike and Bayesian information criterion. Bayesian inference analysis produced a shallow tree without clearly defined groups and with low supports of nodes (Fig. 2). Analyses in AWTY indicated that convergence of the runs was reached. The median-joining network on 41 unique haplotypes revealed a moderately diversified pattern with particular geographic differentiation of particular haplogroups (Fig. 3). Most haplotypes were represented by a single individual. Seven haplotypes were shared by several individuals. First, haplotype CZ-OD5 was shared by four individuals (three of them from identical Krkonoše locality (CZ-OD) one from the disjunct Rhön Mts.) and another six pairs of specimens shared their haplotypes. Distant isolated populations from Rhön Mts. in south Germany and Žd'árské vrchy Mts. in Czech Republic (CZ-KR) cluster with haplotypes from Sudetic Mountains (Krkonoše Mts. and Jizerské hory Mts.), suggesting close relationships among respective areas. Carpathian haplotypes show evident relationship with Sudetic populations (Fig. 3).
Fig. 5.
Bayesian skyline plot showing changes in the effective population size estimate (Ne) overtime measured in mutations per site. The black line depicts the median estimate, and the margins of the blue area represent the highest 95 % posterior density intervals (left). Graph of the mismatch distribution test, which compared the observed frequencies of pairwise differences with the frequencies expected under the demographic model of population expansion (right).
Landscape genetic analysis revealed four groups (Fig. 4). The first group comprises populations from north-western part of study area in the Czech Republic, Austria and Germany (Krkonoše Mts., Šumava Mts., eastern Alps — Landl and isolated Rhön Mts. and Žd'árské vrchy Mts.). Second group (south-western) consist of populations occurring in the south-western part of the study area in Slovenia (south-eastern Alps and northern Dinaric Mts.). Populations from southern Dinaric Mts. in Bosnia and Herzegovina and Montenegro (Bjelašnica Mts., Šator Mts., Zelengora Mts. and Komovi Mts.) belong to the third group (south-western cluster). Fourth group contains populations from the north-east part of the species range (the Carpathians, Fig. 4). The inferred maps of posterior probabilities of population membership are given in supplementary material Fig. S1. Summary statistics of genetic variability were computed for these four groups (Table 2).
Genetic distances among individuals correlated positively with geographic distances. Although isolation by distance was significant, geographic distance explained only tiny portion of genetic heterogeneity (R2 = 0.046, Z = -4 693.332, p < 0.001). Bayesian skyline plot for entire dataset (Fig. 5) as well as for each group (see supplementary material Fig. S2) indicated that Alpine shrew underwent population growth followed by stagnation. Mismatch distribution test retrieved bimodal distribution of pairwise genetic distances (Fig. 5). Statistical evaluation of spatial expansion show that there is no support for stable population, raggedness index = 0.0058 is not significant, (τ = 7.9358, θ = 4.7192, θ1 = 265.9356). Although the results of neutrality tests were negative but not significant, they may point on population growth in the past (Table 2).
Table 2.
Summary statistics of genetic variability for groups resultant from landscape genetic analysis. Abbreviations: standard deviation (SD), significant value (*).
Discussion
Large population in lower altitudes during glacial period
In this study, we provide information about phylogeographic structure of the Alpine shrew, which occurs in the “Alpine archipelago” without long distance disjunctions (Schmitt 2009). We found low differentiation in mitochondrial DNA with evolution of site-specific haplogroups. This pattern contrasts with other soricid species associated with high mountain systems, where deeply divergent genetic lineages were ascertained in habitat islands. For example, allopatric differentiation likely occurred among clades of Sorex bedfordiae species complex, inhabiting high mountain systems in southeastern part of Tibetan Plateau (Chen et al. 2015). Allopatric differentiation of this species is connected with sky island pattern, facilitated by isolation of particular mountain ridges by large river valleys. In contrast, shallow phylogeny of S. alpinus with lack of deep endemic branches implies relatively rapid and recent differentiation from an ancestral population. Similar patterns were determined in several clades of other soricids preferring intermediate elevations as in Sorex monticolus species complex inhabiting the Rocky Mountains and adjacent regions in western North America (Demboski & Cook 2001, Himes & Kenagy 2010) or Sorex cinereus in the Appalachians in eastern North America (Sipe & Browne 2004).
The climatic conditions of glacial periods included both cooling and desiccation of the environment with respect to the inter-glacial periods. Temperature drop resulted in altitudinal shifts of communities. It is likely that during glacial peaks, species that recently avoid habitats above 2500 m a.s.l. occurred ex situ with respect to the current distributional range, that is, in lower altitudes. Drying of the climate caused reduction of the forested area and spread of semi-open and open habitats. Species that exhibit low dispersal abilities and relatively narrow ecological niches likely respond to the climatic oscillations of the Pleistocene by altitudinal shifts rather than long distance movements. For these species, refugia of the coniferous mountain forest, restricted to the foothills of glaciated high mountain systems (Schmitt & Haubrich 2008), may represent suitable conditions. However, existence of multiple centers of glacial survival often results in evolution of two or more deeply separated lineages, which is not the case of S. alpinus. Therefore, we presume that S. alpinus was fairly continuously distributed during the last glacial period. The phylogeographic pattem of shallow divergences among haplotypes corroborates the hypothesis that the glacial range was more extensive than the current one. This perspective is also supported by demographic characteristics. For example, the BSP shows that the population size culminated before recent. However, the bi-modal mismatch distribution could indicate different demographic scenario in different parts of the range. For example, BSP for the Alps shows population growth in recent past, probably concordant with the current interglacial period. This result corroborates our hypothesis about shift of the species to lower altitudes during glacial period. As the Alps were covered with massive glacier during LGM, the ecological conditions of the region were probably not suitable for the species in this time. Nevertheless, evolution of the semifossorial niche in the Alpine shrew is also consistent with the scenario of a more continuous range during the last glacial period, as this species is capable utilizing different depths in the scree habitats which in turn might have helped stenoecious organisms to survive during glacial periods.
Although the distribution of S. alpinus is relatively continuous, gaps within its range exist between the Alps and the Dinarids and between the Sudetes and the Carpathians, with Danube representing the zoogeographical barrier. In both cases, similar genetic lineages inhabit habitat islands on neighboring mountains. This further supports our conclusion that S. alpinus exhibited continuous population during the latest glacial period which was subsequently fragmented and shifted to higher altitudes (Schmitt 2009).
In addition to the high mountain systems of the Alps and Carpathians, the space of Central Europe contains numerous middle high mountains which still host glacial relict populations (Schmitt 2009), including S. alpinus. It is likely that these mountains and adjacent areas with lower altitudes could play substantial role in Vistulian glacial range of S. alpinus. Similar phylogeographical structure to that observed in S. alpinus was ascertained in many other montane organisms. For example, the distribution of the butterfly Erebia epiphron includes the high mountain systems (Alps, Pyrenees) and many isolated range islands in moderate altitudes, including Krkonoše Mts. and Jeseníky Mts. in Central Europe and this species exhibits striking genetic similarity between Jeseníky and some Alpine lineages (Schmitt et al. 2006). Further sampling in southwestern part of the Alpine shrew range will be necessary, however, to fully resolve this issue. It would also be very interesting to examine museum samples from the Pyrenean region where the Alpine shrew is now probably extinct (Mitchell-Jones et al. 1999).
Range breakup in interglacial period and formation of relict enclaves
Moderate diversification and evolution of shallow site-specific genealogy corroborates the hypothesis that the range of S. alpinus became reduced and fragmented during the latest interglacial period. The process of fragmentation of the species suitable habitat was probably gradual. During the Atlantic period of the current interglacial, humid climate was typical for most of the area of Central Europe and contributed to the formation of the closed woodland. Humid environment and presence of closed woodland represent one of the most important habitat requirements of the Alpine shrew (e.g. Kratochvíl & Grulich 1950, Spitzenberger 1978, 2001). Therefore, we believe that suitable habitats for the Alpine shrew existed also outside the extant range. Anthropogenic deforestation of many regions in lower altitudes in Epiatlantic (Horáček & Ložek 1988) might have further facilitated isolation of suitable habitat patches of many glacial relicts, including the Alpine shrew.
Our landscape genetic analyses supported the four cluster model, corresponding with major European high mountain systems. This model fits well with the classical subspecific division of the species, i.e. alpinus in the Alps, hercynicus in Sudetic Mts. and tatricus in the Carpathians (Spitzenberger 1990); populations from the southern Dinarids still lack a subspecific status. Phylogeographic links between high mountain systems and the hilly regions between them were described in many species (Schmitt 2009). The small isolated enclaves of Alpine shrews in discrete hills (Rhön, Žd'árske vrchy) do not show a genetic signal of rapid evolution of distinct lineages. Consequently, there is no evidence for peripatric model. Nevertheless, analysis of fast evolving nuclear loci (i.e. SNPs or microsatellites) will be necessary to describe in more detail the recent population genetic structure of currently isolated populations.
Acknowledgements
This study was supported by the Charles University in Prague, project GA UK No. 575312 and by Czech University of Life Sciences Prague, project IGA 20155014. We are grateful to J. Obuch, I. Baláž, P. Mückstein, M. Švátora, K. Júzová, E. Tošenovský, V. Hlaváč and H. Meinig for providing the tissue samples. We thank T. Ježková (University of Arizona) for improving the English.
Literature
Appendices
Supplementary online material
Fig. S1. Maps of probability of population membership from Geneland analysis, black dots represent sample sites. Subpopulations are displayed by light colours and the contour lines and colour scales denote the values of the posterior probabilities of subpopulation memberships. Particular panels emphasize the following domains: a) Northern Dinarids, b) Southern Dinarids, c) Carpathians, d) Sudetic Mts.
Fig. S2. Bayesian skyline plot for four groups showing changes in the effective population size estimate (Ne) over time measured in mutations per site. The black line depicts the median estimate, and the margins of the blue area represent the highest 95 % posterior density intervals (URL: http://www.ivb.cz/folia_zoologica/supplemetarymaterials/starcova_et_al._fig._s1,_fig._s2.doc).
