The history of Crimean red deer population and Cervus phylogeography in Eurasia

The history of Crimean red deer population and Cervus phylogeography in Eurasia Abstract The present distribution of many species is a result of climatic changes during the Pleistocene and human activity. The impact of climate has been accompanied by restrictions of populations into refugia during glacial periods, and subsequent expansions during more favourable conditions, whereas human influence has been associated with hunting practices and translocations. One mammalian species that has been subject to such transformations is the red deer, Cervus elaphus, but the exact nature of these changes has been difficult to determine using only modern DNA. In this study, we obtained new cytochrome b sequences from subfossil remains of deer found in the Crimean Peninsula. A comparison of these sequences with the available recent and ancient sequences allowed to us to reconstruct phylogeographic relationships between Cervus lineages and to determine their potential migration routes at both local and Eurasian scales. Our analyses showed that the Crimean Peninsula was not a glacial refugium for red deer, but rather that red deer colonized Crimea in three independent waves from both Western and Eastern red deer populations. The immigrations were related to local extinctions and replacements of native populations. ancient DNA, Cervus elaphus, Crimea, phylogeny, phylogeography, Pleistocene, red deer INTRODUCTION During the Late Pleistocene, both environmental change and human activity influenced the distribution and genetic diversity of many species (Hewitt, 1996; Burney & Flannery, 2005; Stewart et al., 2010; Lorenzen et al., 2011; Sandom et al., 2014; Cooper et al., 2015). These factors led to the extinction of many taxa and shaped the modern-day phylogeographic patterns of those that survived. Genetic analyses suggested repeated retreats of ancient populations into refugia during unfavourable conditions and subsequent expansions after a suitable environment stabilized (Hewitt, 1996; Taberlet et al., 1998; Bennett & Provan, 2008). One of the species subjected to such climatic and environmental change was the red deer (Cervus elaphus). Its current distribution is assumed to be strongly influenced by its colonization history during the Late Pleistocene and early Holocene and by human activities throughout the Holocene. Phylogeographic analyses showed two main clades of this species (Kuwayama & Ozawa, 2000; Ludt et al., 2004; Skog et al., 2009). Individuals in Asia and North America have been assigned to the Eastern clade, and European, Middle Eastern and African populations belong to the Western clade. These clades represent populations that are adapted to different climatic conditions (Geist, 1998). The European red deer has a mixed feeding strategy, whereas the Eastern deer is a cold-adapted open-country grazer (Geist, 1998). Estimates of divergence times among the Western and Eastern lineages differ among studies and range from less than a million (Kuwayama & Ozawa, 2000; Mahmut et al., 2002; Polziehn & Strobeck, 2002) to 6–7 Mya (Randi et al., 2001; Ludt et al., 2004). Phylogeographic analyses of European red deer populations have shown that their isolation in Pleistocene refugia had a major impact on their genetic diversity and distribution (Skog et al., 2009; Niedziałkowska et al., 2011). These studies found three mtDNA lineages, marked as haplogroups A, B and C. It has been inferred that these haplogroups correspond to isolation in Iberian, African/Sardinian and Balkan refugia, respectively. Nowadays, haplogroup A is present in Western, Northern and Central Europe. Haplogroup B is restricted to North Africa and Sardinia, whereas haplogroup C is found in Southern and Eastern Europe. Human management over the last centuries is believed to have disrupted the genetic diversity of red deer populations in Europe that arose during the long-term isolation and subsequent post-glacial decolonization. Niedziałkowska et al. (2011) found that haplogroups from lineages A and C occur together in Poland, Belarus and Lithuania. This could be a result of either natural overlap and contact zones or human-mediated translocations that occurred at the end of the 19th and beginning of the 20th centuries. A similar mixture of mtDNA lineages has been found in populations from Bulgaria, Hungary and the Czech Republic (Markov et al., 2012, 2015; Krojerová-Prokešová, Barančeková & Koubek, 2015). Red deer translocations have also been suggested for populations in Germany (Ludt et al., 2004), Scotland (Carden et al., 2012) and Spain (Fernández-García et al., 2013). Moreover, Kuznetsova et al. (2007) found the mtDNA lineages of both Western and Eastern red deer in the Ukrainian population, which probably resulted from the introduction of Cervus elaphus maral to Europe. Therefore, to better understand the impact of climate and environmental conditions on the natural distribution of red deer populations during in the last glaciation, it is crucial to analyse samples from the period before the intensive human management. This can be achieved by ancient DNA (aDNA) analyses, as demonstrated in earlier studies of red deer remains, which indicated that the Late Quaternary history of red deer was more complex than previously assumed (Stankovic et al., 2011; Meiri et al., 2013). Its past genetic diversity was greater than at present, and the distribution of mtDNA lineages was also different. Meiri et al. (2013) noticed that mtDNA lineages only found today in Southern and Eastern Europe were also present in England and Spain before the Last Glacial Maximum (LGM). Analyses of Crimean red deer (Stankovic et al., 2011) indicated that both Western and Eastern lineages were present on this peninsula in the Late Pleistocene, suggesting that conclusions about isolation by distance, which were based on contemporary samples, could be misleading. The Crimean Peninsula has been proposed as a Pleistocene refugium for several species. Biome reconstructions have revealed that forest and steppe-forest vegetation with pine and broadleaved species persisted in Crimea through the last glaciation (Markova, Simakova & Puzachenko, 2009), although Cameron, Pokryszko & Horsák (2013) have argued that full forest cover probably did not survive. Environmental conditions appear to have been stable enough to create a refugium for small mammals, snails and trees (Gerasimenko, 2007; Markova et al., 2009; Markova, 2011; Cameron et al., 2013). Black Sea-level fluctuations in the Pleistocene may periodically have cut off the peninsula from the mainland. The sea level during the Eemian interglacial was few metres higher than today, which could result in an island (Richter, 2005). Crimea became the part of the mainland during colder periods due to sea-level decrease. This level was lower by 30–120 m than present during the Vistulian Interpleniglacial (MIS3), when Crimea was connected with Dobruja, a historical area located in the eastern part of Romania and Bulgaria (Stanko, 2007; Demidenko, 2008). The sea-level fluctuations may therefore have influenced the history of the red deer in the region. Deer populations may have been isolated in the peninsula during warmer events and migrations from adjacent regions may have occurred during colder phases. Similar landscape changes did not happen in other European refugia, which makes Crimea an unusual and interesting area for studying species’ response to Pleistocene climate and geographic changes. In order to test the hypothesis by Stankovic et al. (2011) that the Crimean Peninsula served as Europe’s refugium, from which red deer colonized the Eurasian continent, we obtained additional aDNA from red deer specimens dated to the Late Pleistocene and Holocene. If this hypothesis is correct, we would expect that the Late Pleistocene red deer population would survive the LGM. Therefore, the genetic continuity between samples from before and after the LGM should be observed. Since ancient samples from Crimea represented both Western and Eastern lineages of red deer (Stankovic et al., 2011), we also included all Eurasian available red deer sequences published so far, to reconstruct phylogeographic relationships and determine their potential migration routes at a Eurasian scale. MATERIAL AND METHODS Samples and DNA analyses We analysed 33 samples from red deer remains excavated from ten caves located on three plateaus in the Crimean Mountains (Fig. 1; Supporting Information, Table S1). The majority of the samples were obtained from the Emine-Bair-Khosar Cave (EBK), which is an important palaeontological site rich in vertebrate remains (Vremir & Ridush, 2005; Ridush et al., 2013) (Supporting Information, Fig. S1 and Table S2). Seventeen samples yielding mtDNA were radiocarbon dated in the Poznań Radiocarbon Laboratory and Gliwice Absolute Dating Methods Centre in Poland. Five other samples had been dated previously. The age of the samples ranged from c. 47000 to 900 cal yr BP (Fig. 2; Supporting Information, Table S1), and thus covered a broad time range from the Late Pleistocene to the Holocene. Figure 1. View largeDownload slide A, Crimean Peninsula with the research area (in rectangle) and palaeocoastlines in four time periods. B, research area with the location of caves: 1 – Emine-Bair-Khosar, 2 – K18, 3 – Pavucha, 4 – Primula, 5 – Binbash Koba, 6 – Krapivnyi Grotto, 7 – Dizel, 8 – Bilosnizhka. 9 – Lavrushka, 10 – Sokil. Figure 1. View largeDownload slide A, Crimean Peninsula with the research area (in rectangle) and palaeocoastlines in four time periods. B, research area with the location of caves: 1 – Emine-Bair-Khosar, 2 – K18, 3 – Pavucha, 4 – Primula, 5 – Binbash Koba, 6 – Krapivnyi Grotto, 7 – Dizel, 8 – Bilosnizhka. 9 – Lavrushka, 10 – Sokil. Figure 2. View largeDownload slide Distribution of calibrated radiocarbon dates in time for red deer samples from Crimea. Median (circle) and 95.4% confidence interval (whiskers) are shown. The circles with arrows mean open dates. Phylogenetic clades, to which the samples belong, are shown in parentheses. Figure 2. View largeDownload slide Distribution of calibrated radiocarbon dates in time for red deer samples from Crimea. Median (circle) and 95.4% confidence interval (whiskers) are shown. The circles with arrows mean open dates. Phylogenetic clades, to which the samples belong, are shown in parentheses. The DNA extraction was performed by the Institute of Genetics and Biotechnology of the University of Warsaw, in a laboratory dedicated to aDNA analyses. The amplification of the cytochrome b (cytb) sequences (1140 bp) was performed in multiplex PCR using 12 overlapping primer pairs described in Stankovic et al. (2011). Full details regarding DNA extraction and amplification are given in Supporting Information. Multiplex PCR products were used to prepare libraries for pyrosequencing following the protocol from Stiller et al. (2009) and were sequenced on a Roche 454 sequencing platform at the Institute of Biochemistry and Biophysics, PAS. The resulting sequences were trimmed in Bioedit v. 7.2.5 (Hall, 1999) and assembled into contigs using SeqMan Pro software (DNASTAR, Inc.). To ensure authenticity of the sequence data, each sample was amplified in two independent multiplex PCRs and contigs from both replicates were used to create consensus sequences according to guidelines proposed by Stiller et al. (2009). Phylogenetic analyses Data sets were obtained by comprehensive searches of the GenBank database for cytochrome b homologs in members of the Cervini tribe. Many of the red deer sequences studied here, including fossil samples, were obtained by various research groups, such as Ludt et al. (2004), Skog et al. (2009), Kuznetsova, Danilkin & Kholodova (2012) and Meiri et al. (2013, 2014). Phylogenetic analyses were based on two alignments comprising non-redundant sequences: a longer set of 190 sequences, 1140 bp in length and a shorter set of 175 sequences with 423 bp. We applied five approaches for phylogenetic reconstruction: Bayesian inference in MrBayes (Ronquist et al., 2012) and PhyloBayes (Lartillot, Lepage & Blanquart, 2009), maximum likelihood analyses in TreeFinder (Jobb, von Haeseler & Strimmer, 2004) and morePhyML (Criscuolo, 2011) based on PhyML (Guindon et al., 2010), and the neighbor-joining method in PAUP* (Swofford, 1998). In the MrBayes analyses, we assumed separate mixed + I + Γ(5) models for each of three codon positions and two independent runs, each using 32 and 16 Markov chains for the longer and shorter sets, respectively. Trees were sampled every 100 generations for 10000000 generations. In the final analysis, we selected trees from the last 3062000 (the longer set) and 4664000 (the shorter set) generations that reached convergence. In the PhyloBayes approach, two independent Markov chains were run through 100000 (the longer set) and 200000 (the shorter set) cycles assuming the CAT GTR + Γ(5) model. After obtaining convergence, the last 50000 and 10000 trees from each chain were collected to compute the posterior consensus, respectively. In TreeFinder, we applied a search depth set to 2, separate substitution models for each of three codon positions according to the TreeFinder Propose Model (Supporting Information, Table S3). Trees inferred with morePhyML and PAUP were based on the best-fit substitution models found in jModeltest (Darriba et al., 2012) among 1624 candidate models (Supporting Information, Table S3). The best search algorithm NNI + SPR was applied in morePhyML. To assess the significance of branches, non-parametric bootstrap analyses were performed on 1000 replicates in TreeFinder, PhyML and PAUP. Molecular dating Divergence times were estimated with the software BEAST (Drummond et al., 2012) on the full-length alignment comprising 216 sequences, using as many samples as possible with assigned dates and geographic localizations. We used a fixed MrBayes tree that was inferred according to the procedure described above, but with 32 Markov chains and the last 3861000 generations to summarize the final tree. Based on that, we created an ultrametric tree using Sanderson’s approach (Sanderson, 2002) in the chronopl function in R (R Core Team, 2014). We applied separate substitution models for each of three codon positions selected in jModeltest among 1624 candidate models (Supporting Information, Table S2). The age of the root of the tree was set as a uniform prior with the lower and upper limits 5.1 and 7.3 Mya. This age corresponds to the separation time of Axis and Rucervus from other Cervini based on three methods (Hassanin et al., 2012). Moreover, we used 25 tip dates associated with sequences obtained from fossil and historical samples. To unify the dates, we applied the same calibration procedure in OxCal 4.2.3 (Ramsey, Scott & van der Plicht, 2013) using curve IntCal 13 (Reimer et al., 2013) for original uncalibrated dates. We tested strict and uncorrelated lognormal relaxed clock models with a coalescent constant size tree and separate clock models for particular data partitions. Finally, we selected the relaxed clock model because the SD of the relaxed clock revealed significant variation among branches. Posterior distributions of parameters were estimated for 50000000 generations with a sampling frequency of 1000 steps. The convergence and sufficient sampling were checked using Tracer (Rambaut et al., 2014). The phylogenetic trees were summarized in TreeAnnotator (Drummond et al., 2012) with a 10% burn-in. Approximate Bayesian computation Based on archaeological and historical data, we defined a set of the 12 most likely models that described the temporal pattern in genetic diversity of Crimean red deer (Fig. 3). The best-fitted model was selected by approximate Bayesian computation (ABC) (Beaumont, Zhang & Balding, 2002). For defining the models, we considered the following scenarios: Figure 3. View largeDownload slide Population models considered in approximate Bayesian computation analyses. The most likely model H assuming 5.14 and 10% mutation rate per Myr was outlined by the red rectangle. Figure 3. View largeDownload slide Population models considered in approximate Bayesian computation analyses. The most likely model H assuming 5.14 and 10% mutation rate per Myr was outlined by the red rectangle. Two possible scenarios regarding the pre-LGM origin of Crimean population: one including a long-standing panmictic population (Fig. 3A, C, E, G, I, K) and the second assuming an admixture of Western and Eastern populations (Fig. 3B, D, F, H, J, L); Two possible scenarios regarding the post-LGM (until c. 2000 cal yr BP) origin of the Crimean population: in one scenario, Crimea was populated by deer that survived locally (Fig. 3A–F). In the other, Crimea was colonized by Western European lineages (Fig. 3H–L); Three scenarios regarding the present origin of the Crimean population: in one scenario, the present population was derived from the population inhabiting Crimea before 2000 cal yr BP (Fig. 3A, B, K, L); in the second scenario, the present population split from a common ancestor with the population that colonized Crimea after the LGM and until 2000 cal yr BP. Here we considered two split times: before the LGM (Fig. 3I, J) and after the LGM (Fig. 3C, D). In the third scenario, the present population was colonized from Western European lineages that were unrelated to the post-LGM (till c. 2000 cal yr BP) population (Fig. 3E–H). The 12 tested models resulted from combining these scenarios. Simulations were performed under a coalescent framework, where unknown parameters and their uncertainness were assessed by Bayesian prior probability distributions. We carried out two rounds of pilot simulations to optimize the parameters. The priors for effective population size (Ne) had an exponential distribution in pilot simulations to better screen an ample range. Then, we replaced them with uniform distributions in the final simulations. The times to population splits or demographic changes were set by uniform priors as well. Uncertainty in radiocarbon dating of sample ages was included by assigning normal priors (except in four cases where exponential priors were used because of infinite radiocarbon dates; see Supporting Information). We employed a 4-year generation time, based on the literature (Pérez-Espona et al., 2009; Rosvold et al., 2012) and a transition/transversion (Tr/Ts) bias of 0.945, a gamma shape parameter of 0.05 and a mutation rate of 5.14% per Myr (Skog et al., 2009). We calculated the basic parameters: the number of haplotypes, private haplotypes and segregating sites and nucleotide diversity for three sets of samples grouped by time periods: pre-LGM, post-LGM to 2000 cal yr BP and the present. In addition, we employed the number of pairwise differences and FST among pairs of consecutive samples (Supporting Information, Table S4). We ran 2000000 simulations for each model in the final analyses. The best-fitting model was selected with Bayes factors (BFs) that were calculated as ratios of model likelihoods between pairs of models for all 66 pairwise comparisons of the 12 models (Supporting Information, Table S5). The model likelihoods were estimated by (1) a direct approach where they were directly equalled to their relative acceptance ratios, (2) a weighted approach, which also employs the acceptance ratios after being weighted with an Epanechnikov kernel evaluated in the Euclidean distances to the observed data and (3) a logistic regression in which summary statistics represent explanatory variables of model likelihoods (Fagundes et al., 2007). The consistency of the model likelihoods and the BFs was assessed by applying the procedure with 10 different acceptance proportions (from 0.001 to 0.01%). Since the mutation rate, and possible dependency on temporal scale (Ho, 2014), could impact our estimated parameters and model choice, we tested the robustness of our results with a range of different mutation rates: 10, 20 and 30% per Myr. Finally, to assess the probability of correct model selection, we estimated the statistical power of the analyses with pseudo observed data sets (PODs), where simulated data sets play the role of real ones. We generated 1000 PODs, which were analysed with the same reference tables and 2000000 simulations each. Simulations, model choice and POD analyses were performed in the software BaySICS (Sandoval-Castellanos, Palkopoulou & Dalén, 2014), whereas the Tr/Ts bias and gamma parameter were calculated in MEGA 6 (Tamura et al., 2013). RESULTS Phylogenetic relationships within red deer lineages including Crimean deer sequences We successfully amplified and obtained cytochrome b sequence from 24 samples, which gave us a 73% success rate of DNA extraction for the subfossil samples. We identified 16 haplotypes. None of the negative controls yielded DNA. The sequences were aligned to all available cytochrome b sequences from Cervini members including ancient and modern samples of C. elaphus. The five applied approaches produced similar trees. Deer sequences assigned to the Cervus genus do not form a monophyletic clade but are split into Western and Eastern subclades, with strong and moderate support, respectively (Supporting Information, Figs S2, S3). The Eastern subclade is very significantly grouped with the Thorold’s deer Przewalskium albirostris sequences. All newly obtained sequences from the Crimean samples group significantly either with the Western or Eastern lineages (detailed phylogenies in Figs 4, 5; Supporting Information, Figs S4, S5). To study the detailed relationships between Cervus sequences, we considered the shorter set (Supporting Information, Figs S4, S5) in addition to the longer one (Figs 4, 5) because it contained important subfossil samples. Figure 4. View largeDownload slide MrBayes tree based on the longer alignment for the Eastern lineage of red deer (the clade E in Supporting information, Fig. S2). Clades E2 and E3 represent C. elaphus canadensis (wapiti), whereas E4–E7 C. elaphus nippon (sika). Sequences newly obtained in this study from Crimea specimens are in bold. The asterisk indicates that a given sequence is represented by more than one accession number. The shaded circles at nodes correspond to the number of methods (out of five) that recovered the given node. Probable age, time period or calibrated dates of samples are shown in square brackets. Numbers at nodes, in the order shown, correspond to posterior probabilities estimated in MrBayes (MB) and PhyloBayes (PB), bootstrap values obtained in PhyML (PH), TreeFinder (TF) and PAUP by the neighbor-joining (NJ) method. Values of the posterior probabilities and bootstrap percentages ≤ 0.50 and 50%, respectively, were omitted or indicated by a dash ‘–’. Figure 4. View largeDownload slide MrBayes tree based on the longer alignment for the Eastern lineage of red deer (the clade E in Supporting information, Fig. S2). Clades E2 and E3 represent C. elaphus canadensis (wapiti), whereas E4–E7 C. elaphus nippon (sika). Sequences newly obtained in this study from Crimea specimens are in bold. The asterisk indicates that a given sequence is represented by more than one accession number. The shaded circles at nodes correspond to the number of methods (out of five) that recovered the given node. Probable age, time period or calibrated dates of samples are shown in square brackets. Numbers at nodes, in the order shown, correspond to posterior probabilities estimated in MrBayes (MB) and PhyloBayes (PB), bootstrap values obtained in PhyML (PH), TreeFinder (TF) and PAUP by the neighbor-joining (NJ) method. Values of the posterior probabilities and bootstrap percentages ≤ 0.50 and 50%, respectively, were omitted or indicated by a dash ‘–’. Figure 5. View largeDownload slide MrBayes tree based on the longer alignment for the Western lineage of red deer (the clade W in Supporting Information, Fig. S2). Sequences newly obtained in this study from Crimean specimens are shown in bold. For further explanation, see Figure 4. Figure 5. View largeDownload slide MrBayes tree based on the longer alignment for the Western lineage of red deer (the clade W in Supporting Information, Fig. S2). Sequences newly obtained in this study from Crimean specimens are shown in bold. For further explanation, see Figure 4. In the Eastern lineage, two main well-supported groups representing wapiti, Cervus canadensis, and sika deer, Cervus nippon, can be recognized (Fig. 4; Supporting Information, Fig. S4). The wapiti group consists of clades marked as E1 (only in the tree for the smaller set, Supporting Information, Fig. S4), E2 and E3. The basal clade E1 comprises exclusively ancient sequences from Alaska and north-eastern Siberia (Meiri et al., 2014). The pre-LGM samples from north-eastern Siberia (Meiri et al., 2014) are located in the E2 and E3 clades too. The E2 clade also includes other samples from Central and central-eastern Asia. The deer samples in the E3 clade show a much wider distribution. Besides Central Asia, north-eastern Siberia and South Korea, they were also found in North America and Crimea, where they are represented by only pre-LGM samples. In the tree for the longer set, all Crimean sequences branch off at the base of E3 clade, but for the shorter set, they are placed inside the clade without any support (Fig. 4; Supporting Information, Fig. S4). In the sika deer group, four clades can be recognized (Fig. 4; Supporting Information, Fig. S4) including Asian (clades E4 and E5) and Japanese deer separated into its ‘northern’ (clade E6) and ‘southern’ (clade E7) (Nagata et al., 1999; Groves, 2005). Deer introduced to the Czech Republic are also present in these clades. However, no Crimean samples were reported among the sika deer clade. The Western red deer lineage is diversified into nine distinct clades, which are inferred by at least three out of the five methods applied (Fig. 5; Supporting Information, Fig. S5). Sequences from Crimean deer are distributed among four clades (W5, W6, W7 and W9). The Late Pleistocene samples from Crimea (post-LGM) are present in clade W5, which also contains other sequences from the southern part of Eastern Europe: modern samples from Krasnodarsky Krai, Voronezh and the Caucasus (from the pre-LGM period) (Meiri et al., 2013). Much older Crimea sequences from the pre-LGM are grouped with sequences from the Balkans, Italy, Central Europe and Turkey in the W6 clade (corresponding to haplogroup C). In the smaller set tree (Supporting Information, Fig. S5), this clade also includes pre-LGM specimens from the Caucasus and the Urals (Meiri et al., 2013) as well as a modern Crimean sequence. In the longer alignment tree (Fig. 5), the W7 clade contains exclusively post-LGM Crimean sequences, whereas in the shorter set tree (Supporting Information, Fig. S5), it includes additionally a modern sequence from Italy, two post-LGM sequences from Serbia and two ancient sequences from England, with pre- and post-LGM ages (Meiri et al., 2013), respectively. The W9 clade (corresponding to haplogroup A) has the highest sequence diversity. Besides Mediaeval (eighth to 11th centuries) and modern sequences from Crimea, it comprises deer mainly from Western, Central and Northern Europe including Great Britain, and sequences from regions far to the east, in Russian Karelia and Belarus, and deer introduced to South Korea and New Zealand. In the tree based on the short alignment (Supporting Information, Fig. S5), the basal position of this clade is occupied by a pre-LGM sample from Belgium (Meiri et al., 2013). Post-LGM specimens from Spain are branched off next. The deer from Crimea was not found in the basal clades of the Western lineage: W1 including Cervus elaphus yarkandensis and bactrianus from China and Uzbekistan; W2 (in the shorter alignment tree only, Supporting Information, Fig. S5) containing exclusively post-LGM sequences from Spain (Meiri et al., 2013); W3 (in the shorter alignment tree only, Supporting Information, Fig. S5) with a historical Jordan sample (Meiri et al., 2013); and W4 comprising the sequences of C. elaphus maral from Turkey and Iran. The W8 clade (corresponding to haplogroup B) is also absent from Crimean samples and includes deer from Tunisia, Sardinia, Hungary, the Swiss Alps, a deer introduced to Germany and pre-LGM samples from the Caucasus – Supporting Information, Figure S5. Population history of Crimean deer In the ABC model choice analysis, the same model H was selected for mutation rates of 5.14 and 10% per Myr (Fig. 3H; Table 1; Supporting Information, Table S5). This model assumed that the pre-LGM population was an admixture of Eastern (Asian) and Western (Europe) lineages, whereas the two post-LGM populations (post-LGM to 2000 yr BP, and present) resulted from two unrelated waves of colonization that originated in the West before 50 ka. The BF of the chosen model was relatively large (Table 1; Supporting Information, Table S5). For mutation rates of 20 and 30% per Myr, model G instead obtained the highest support (Supporting Information, Tables S6 and S7). The model G assumed that the pre-LGM Crimean population was panmictic and long standing. Table 1. Model likelihoods and Bayes factors among the 12 models compared by approximate Bayesian computation Model Model likelihoods A B C D E F G H I J K L A 0.0006 – 0.348 0.071 0.128 0.021 0.021 0.004 0.001 0.006 0.003 0.095 0.034 B 0.0018 2.871 – 0.205 0.368 0.059 0.059 0.012 0.004 0.017 0.009 0.272 0.097 C 0.0087 14.000 4.876 – 1.793 0.289 0.288 0.056 0.020 0.085 0.044 1.327 0.474 D 0.0049 7.806 2.719 0.558 – 0.161 0.161 0.031 0.011 0.048 0.024 0.740 0.264 E 0.0300 48.419 16.865 3.459 6.202 – 0.996 0.194 0.068 0.295 0.151 4.590 1.640 F 0.0302 48.613 16.933 3.472 6.227 1.004 – 0.195 0.068 0.296 0.152 4.609 1.647 G 0.1547 249.355 86.854 17.811 31.942 5.150 5.129 – 0.348 1.519 0.779 23.639 8.448 H 0.4439 715.903 249.360 51.136 91.707 14.785 14.727 2.871 – 4.362 2.237 67.869 24.255 I 0.1018 164.129 57.169 11.724 21.025 3.390 3.376 0.658 0.229 – 0.513 15.560 5.561 J 0.1985 320.032 111.472 22.859 40.996 6.610 6.583 1.283 0.447 1.950 – 30.339 10.843 K 0.0066 10.548 3.674 0.753 1.351 0.218 0.217 0.042 0.015 0.064 0.033 – 0.357 L 0.0183 29.516 10.281 2.108 3.781 0.610 0.607 0.118 0.041 0.180 0.092 2.798 – Model Model likelihoods A B C D E F G H I J K L A 0.0006 – 0.348 0.071 0.128 0.021 0.021 0.004 0.001 0.006 0.003 0.095 0.034 B 0.0018 2.871 – 0.205 0.368 0.059 0.059 0.012 0.004 0.017 0.009 0.272 0.097 C 0.0087 14.000 4.876 – 1.793 0.289 0.288 0.056 0.020 0.085 0.044 1.327 0.474 D 0.0049 7.806 2.719 0.558 – 0.161 0.161 0.031 0.011 0.048 0.024 0.740 0.264 E 0.0300 48.419 16.865 3.459 6.202 – 0.996 0.194 0.068 0.295 0.151 4.590 1.640 F 0.0302 48.613 16.933 3.472 6.227 1.004 – 0.195 0.068 0.296 0.152 4.609 1.647 G 0.1547 249.355 86.854 17.811 31.942 5.150 5.129 – 0.348 1.519 0.779 23.639 8.448 H 0.4439 715.903 249.360 51.136 91.707 14.785 14.727 2.871 – 4.362 2.237 67.869 24.255 I 0.1018 164.129 57.169 11.724 21.025 3.390 3.376 0.658 0.229 – 0.513 15.560 5.561 J 0.1985 320.032 111.472 22.859 40.996 6.610 6.583 1.283 0.447 1.950 – 30.339 10.843 K 0.0066 10.548 3.674 0.753 1.351 0.218 0.217 0.042 0.015 0.064 0.033 – 0.357 L 0.0183 29.516 10.281 2.108 3.781 0.610 0.607 0.118 0.041 0.180 0.092 2.798 – The Bayes factors have a row to column interpretation (the value at row i and column j is BFij = MLi/MLj). These values correspond to the analysis with a mutation rate of 5.14% per Myr, after averaging the five smallest acceptance proportions (0.001–0.005%) of the closest simulated data sets to the observed one. The best-fit model is shown in bold. View Large Table 1. Model likelihoods and Bayes factors among the 12 models compared by approximate Bayesian computation Model Model likelihoods A B C D E F G H I J K L A 0.0006 – 0.348 0.071 0.128 0.021 0.021 0.004 0.001 0.006 0.003 0.095 0.034 B 0.0018 2.871 – 0.205 0.368 0.059 0.059 0.012 0.004 0.017 0.009 0.272 0.097 C 0.0087 14.000 4.876 – 1.793 0.289 0.288 0.056 0.020 0.085 0.044 1.327 0.474 D 0.0049 7.806 2.719 0.558 – 0.161 0.161 0.031 0.011 0.048 0.024 0.740 0.264 E 0.0300 48.419 16.865 3.459 6.202 – 0.996 0.194 0.068 0.295 0.151 4.590 1.640 F 0.0302 48.613 16.933 3.472 6.227 1.004 – 0.195 0.068 0.296 0.152 4.609 1.647 G 0.1547 249.355 86.854 17.811 31.942 5.150 5.129 – 0.348 1.519 0.779 23.639 8.448 H 0.4439 715.903 249.360 51.136 91.707 14.785 14.727 2.871 – 4.362 2.237 67.869 24.255 I 0.1018 164.129 57.169 11.724 21.025 3.390 3.376 0.658 0.229 – 0.513 15.560 5.561 J 0.1985 320.032 111.472 22.859 40.996 6.610 6.583 1.283 0.447 1.950 – 30.339 10.843 K 0.0066 10.548 3.674 0.753 1.351 0.218 0.217 0.042 0.015 0.064 0.033 – 0.357 L 0.0183 29.516 10.281 2.108 3.781 0.610 0.607 0.118 0.041 0.180 0.092 2.798 – Model Model likelihoods A B C D E F G H I J K L A 0.0006 – 0.348 0.071 0.128 0.021 0.021 0.004 0.001 0.006 0.003 0.095 0.034 B 0.0018 2.871 – 0.205 0.368 0.059 0.059 0.012 0.004 0.017 0.009 0.272 0.097 C 0.0087 14.000 4.876 – 1.793 0.289 0.288 0.056 0.020 0.085 0.044 1.327 0.474 D 0.0049 7.806 2.719 0.558 – 0.161 0.161 0.031 0.011 0.048 0.024 0.740 0.264 E 0.0300 48.419 16.865 3.459 6.202 – 0.996 0.194 0.068 0.295 0.151 4.590 1.640 F 0.0302 48.613 16.933 3.472 6.227 1.004 – 0.195 0.068 0.296 0.152 4.609 1.647 G 0.1547 249.355 86.854 17.811 31.942 5.150 5.129 – 0.348 1.519 0.779 23.639 8.448 H 0.4439 715.903 249.360 51.136 91.707 14.785 14.727 2.871 – 4.362 2.237 67.869 24.255 I 0.1018 164.129 57.169 11.724 21.025 3.390 3.376 0.658 0.229 – 0.513 15.560 5.561 J 0.1985 320.032 111.472 22.859 40.996 6.610 6.583 1.283 0.447 1.950 – 30.339 10.843 K 0.0066 10.548 3.674 0.753 1.351 0.218 0.217 0.042 0.015 0.064 0.033 – 0.357 L 0.0183 29.516 10.281 2.108 3.781 0.610 0.607 0.118 0.041 0.180 0.092 2.798 – The Bayes factors have a row to column interpretation (the value at row i and column j is BFij = MLi/MLj). These values correspond to the analysis with a mutation rate of 5.14% per Myr, after averaging the five smallest acceptance proportions (0.001–0.005%) of the closest simulated data sets to the observed one. The best-fit model is shown in bold. View Large Four models obtained higher likelihood support than the rest (Figs 2I, 3G–J). In all four models, the most recent populations originated before the LGM in two waves of colonization. These models had BF values in the range of 7.5–360, 4.3–3800, 8.6–1560 and 18.0–105 for mutation rates of 5.14, 10, 20 and 30% per Myr, respectively. These results provide strong support for the presence of at least three waves of sequential colonization (and extinction/replacement) of Crimean red deer by unrelated Western populations, and reject that the present-day population originated from a local pre-LGM population or even a post-LGM population that existed until 2000 yr BP. For mutation rates of 5.14 and 10% per Myr, the BF was > 2.3 for the models with the admixed pre-LGM population (Fig. 3H, J), and it was > 1.5 for the models G and I where the pre-LGM population was panmictic (Fig. 3G, I). However, for mutation rates of 20 and 30% per Myr, the converse occurred (BF = 0.12 and 0.06, respectively). Comparison of models also indicates that the support for a Western origin of post-LGM populations in at least two colonization waves is strong and consistent. Finally, the POD analysis showed low statistical power for several models (Supporting Information, Table S8) ranging from 0.22 to 0.48, which rose to 0.28–0.83 when only BFs equal to or larger than the observed ones were taken into account. The lack of statistical power resulted from ‘competition’ between similar models. An exploratory analysis showed that if similar models (e.g. models of Fig. 3G, H) were mixed into one, the statistical power rose to over 0.9 even for the least supported scenarios. Therefore, the POD analysis showed that the model choice procedure is able to separate correctly between scenarios, although it is less sensitive to distinguishing between very similar models. DISCUSSION Phylogeography of Cervus lineages including Crimean samples We analysed the most sequence-rich set of cytochrome b from different Cervini using the greatest variety of methods so far. The obtained global tree topologies are generally in agreement with results of other authors (Pitra et al., 2004; Hassanin et al., 2012) – Supporting Information, Figure S2. Divergence times for particular lineages of Cervini estimated by us are close to those by Gilbert, Ropiquet & Hassanin (2006), although slightly younger (Supporting Information, Fig. S3). Our phylogenetic study clearly recovered two monophyletic clades of Cervus, corresponding to the eastern and western part of its distribution (Randi et al., 2001; Ludt et al., 2004; Pitra et al., 2004; Skog et al., 2009). The analyses of a large number of deer samples, including fossils from Crimea, also enabled us to infer branching order and migration routes of Cervus in more detail. The significant grouping of the Eastern lineage with Przewalskium from China suggests that this lineage emerged somewhere in Asia. The phylogenies are based on an mtDNA marker, which reflects the history of the maternal lineage. Therefore, nuclear markers should be studied to verify the inferred relationships. About 1 Mya, the Eastern lineage separated into sika (C. nippon) and wapiti (C. canadensis) clades, which further differentiated in Central and Eastern Asia (Fig. 4; Supporting Information, Figs S3, S4). The former colonized Japanese islands (Nagata et al., 1999; Groves, 2005) and the latter reached north-eastern Siberia in the pre-LGM period (Meiri et al., 2014). Wapiti deer carrying mtDNA haplotypes belonging to the E3 clade showed the most widespread distribution because they not only colonized North America but also migrated as far west as Crimea, where they survived to ~35 ka. This implies long-distance migrations of the eastern deer, originally confined into Central and Eastern Asia. The colonization of North America was a single event (Meiri et al., 2014) (Fig. 4) and probably started from north-eastern Siberia (Supporting Information, Fig. S3). Our molecular dating indicates that the migration occurred between 12.4 and 75.9 ka. The deer that crossed the land bridge across the present Bering Strait to North America carried at least two mitochondrial haplotypes, represented by clades E1 and E3. The ancestor of the Western lineage most likely lived in Central Asia, as indicated by the earliest diverged W1 clade (Fig. 5; Supporting Information, Fig. S5). This clade includes Yarkand deer (C. elaphus yarkandensis) and Bactrian deer (C. elaphus bactrianus) from China and Uzbekistan, closely related to C. elaphus hanglu from Indian Kashmir (Lorenzini & Garofalo, 2015). Lorenzini & Garofalo (2015) proposed to raise these three deer to species level as Cervus hanglu. The migration of red deer to Europe happened likely about 1 Mya, as implied by the oldest fossils of C. elaphus (Brugal & Croitor, 2007) – Figure 6A. This is supported by our phylogenies, where an early separated clade with Spanish samples (W2) is basal to the rest of the sequences (Supporting Information, Fig. S5). This is also consistent with our estimated split time between clade W1 and the rest: 1.08 (0.57–1.72) Mya (Supporting Information, Fig. S3). The ancient Spanish samples are dated to post-LGM and do not group with other Spanish sequences, which separated much later (the W9 clade), suggesting that they represent an old lineage that became extinct without descendants. A similar placement of these sequences was obtained in the phylogenetic tree by Meiri et al. (2013). The migration route to Europe probably extended across the Near East and Asia Minor (Fig. 6A) because sequences from Iran, Jordan and Turkey clades are closely related with those from Spain (Supporting Information, Fig. S5). Figure 6. View largeDownload slide Inferred migration routes of red deer in Europe. A, the first settlement of Europe from Central Asia across Near East and Asia Minor as well as migrations to the southern region of Eastern Europe including the Crimean Peninsula. B, further migrations from the latter region to Europe and towards the Ural Mountains. C, post-LGM recolonization of the rest territory of Europe including the Crimean Peninsula from Iberian refugium. The role of the Balkan and Apennine refugia is not shown. Figure 6. View largeDownload slide Inferred migration routes of red deer in Europe. A, the first settlement of Europe from Central Asia across Near East and Asia Minor as well as migrations to the southern region of Eastern Europe including the Crimean Peninsula. B, further migrations from the latter region to Europe and towards the Ural Mountains. C, post-LGM recolonization of the rest territory of Europe including the Crimean Peninsula from Iberian refugium. The role of the Balkan and Apennine refugia is not shown. The next diverged W5 clade includes sequences from the Caucasus, Krasnodarsky Krai, Voronezh and Crimea, which suggests that populations from the Near East also migrated north along the western or/and eastern seashore of the Caucasus and inhabited the southern part of Eastern Europe, including Crimea (Fig. 6A). The first migration may have occurred around 0.6 Mya because the oldest red deer from the Caucasus is dated to 583 ka (Blackwell et al., 2005). This agrees with our calculation of separation time between the lineage including Turkey and Iran sequences and other clades, that is 512 (265–832) ka (Supporting Information, Fig. S3). Samples from Crimea and the Caucasus are separated into different clades with short internal branches, which may indicate a quite diversified ancestral population or several rapid migration events to the north. Some deer from the region migrated to the northeast (Fig. 6A) as indicated by ancient sequences from the Urals, identical to those from the Caucasus (Supporting Information, Fig. S5). At least four westward migrations could have occurred from the southern part of Eastern Europe (Fig. 6B) or different haplotypes were carried by a larger population during a smaller number of migration events. This is suggested by clades W6–8 including pre-LGM sequences from Crimea and the Caucasus and from different regions of Europe (Fig. 5; Supporting Information, Fig. S5). In the W6 and W8 clades, these sequences occupy a basal position to others (Supporting Information, Fig. S5). The separation of these clades appears to have begun roughly 468 (244–758) ka (Supporting Information, Fig. S3). Drastic cooling during the LGM (23–19 ka) caused the distribution of red deer in Europe to be limited to the Balkan, Italian and Iberian Peninsulas, south-western France and the region east of the Carpathians (Banks et al., 2008; Sommer et al., 2008; Skog et al., 2009; Sommer & Zachos, 2009; Meiri et al., 2013). After the LGM with the start of Bølling-Allerød warming (15–13 ka), deer expanded from refugial areas into other regions of Europe. Our results are in agreement with those studies. In clade W9 (corresponding to haplogroup A), the basal position is occupied by Spanish sequences dated from c. 19 ka to the Neolithic (Meiri et al., 2013) – Supporting Information, Figure S5. Other Spanish sequences are younger and are intermingled in the tree with sequences from Western and Central Europe, in line with expectations under a scenario where the red deer that colonized other regions of Europe originated in the Iberian Peninsula (Ludt et al., 2004; Sommer et al., 2008; Skog et al., 2009; Sommer & Zachos, 2009; Niedziałkowska et al., 2011; Meiri et al., 2013) – Figure 6C. Sequences from Crimean deer are separated in the W9 clade, which suggests two independent waves of migration into the Crimean Peninsula or one colonization event by a deer population that carried different haplotypes. A subcluster in clade W6 (corresponding to C haplogroup) with sequences from Southern and Central Europe may correspond to the Balkan refugium. Sequences from Turkey are also present in this clade, which could imply a migration between Europe and Asia Minor. Alternatively, deer from the Balkan refugium could have spread further to the south-east. The W8 clade (corresponding to B haplogroup), besides harbouring sequences from Corsican and Tunisian red deer (C. elaphus corsicanus and C. elaphus barbarus), also includes sequences from Hungary and the Swiss Alps. These most likely represent human translocations from Italy. Detailed population history of Crimean red deer The goal of our study was to investigate the role of the Crimean Peninsula as a possible Eastern European refugium. In agreement with Stankovic et al. (2011), we found that red deer from both the Eastern and Western lineages were present in the Crimean Peninsula before the LGM, even though these populations have otherwise been isolated since their divergence (Kuwayama & Ozawa, 2000; Mahmut et al., 2002). Our ABC analyses, assuming the mutation rates of 5.14 and 10% per Myr, showed that the pre-LGM Crimean population was probably the result of a recent mixture of these lineages. However, for higher mutation rates, the best-supported model (G) assumes that the pre-LGM population was panmictic. This discrepancy results from the model selection procedure, where the POD analysis had difficulties distinguishing between very similar models. Based solely on the ABC results, we cannot therefore be certain which of the scenario regarding the pre-LGM population is more probable. However, the scenario presented in model H is in agreement with other results: the early divergence estimate for the two red deer lineages (c. 1.54 Mya), the distribution of pre-LGM specimens in separate clades, E3 and W6 from the Eastern and Western lineages, and radiocarbon dates and the distribution of red deer samples in cave layers. The specimens from clade W6 were usually found in the lower part of sediments in Emine-Bair-Khosar Cave (Supporting Information, Fig. S1). This layer (I) was dated with pollen to the Eemian stage (MIS5e) when the climate was warmer than now (Avdeienko, 2015) – Supporting Information, Table S2. The oldest deer remains (J46 and J47) with 14C radiocarbon dates close or above the upper limit of this method also belong to clade W6. On the other hand, samples from the eastern population, located in clade E3, were frequent in the upper layer (H), dated to the Middle Pleniglacial (MIS3) by 14C and pollen (Avdeienko, 2015) – Supporting Information, Figure S1 and Table S2. The deer lived in a boreal climate, which was cooler than now, and became extinct before the LGM when the accumulation of allochthonous material in the cave stopped. The youngest pre-LGM sample J28 (28332 cal yr BP) again belongs to clade W6. The data suggest that the deer from clade W6 dominated in Crimea earlier, next the deer from clade E3 appeared and again the deer from clade W6 returned (Fig. 2). However, some radiocarbon dates of samples from these two lineages overlap. Therefore, we cannot exclude that the different deer could interact in some periods (e.g. between 47400 and 43700 cal yr BP), although there were also periods when only one deer clade existed, for example before c. 50000 cal yr BP or between 37300 and 34900 cal yr BP. We did not observe a genetic continuity between pre- and post-LGM samples from Crimea (Fig. 2), which suggests a local extinction hypothesis. This is also supported by the lack of red deer remains from the LGM period in Emine-Bair-Khosar Cave, where such remains are generally abundant (Vremir & Ridush, 2005). However, it cannot be excluded that bone accumulation was impossible in the MIS 2 due to a snow-ice plug that covered the entrance to the cave (Ridush et al., 2013). In this time, the Scandinavian Ice Sheet reached its maximum extent, the climate became colder and drier, the sea level dropped and many temperate species retreated to refugial areas in the southern peninsulas (Hewitt, 2004; Hubberten, 2004; Clark et al., 2009). Environmental conditions in Crimea also deteriorated at the end of MIS 3. Pollen, palaeontological and stable isotope analyses indicated that the climate in the LGM was cold and dry, which reduced broadleaved trees and extended steppes (Gerasimenko, 2007; Markova, 2011; Gąsiorowski et al., 2014). Those unfavourable conditions could have contributed to the extinction of the pre-LGM red deer population. However, the results of Markova (2011) revealed that despite climate deterioration, Crimea was a refugium for small mammals. Therefore, we cannot exclude that climate was the only factor that influenced the red deer extinction. The post-LGM samples are distributed into three clades that belong to the Western lineage and were not descendants of the pre-LGM population, as is indicated by our phylogenetic and ABC analyses. Thus, the pre-LGM population became extinct and was replaced by immigrants. Our data suggest that the population replacement started around 27–14 ka when the pre-LGM deer disappeared from Crimea (Fig. 2). The oldest immigrant population is represented by samples J16 and J27 (c. 13900 cal yr BP), which were probably derived from Russian and Caucasian populations (assigned in Fig. 5 as J16* and in Supporting Information, Fig. S5 as JX966135* in clade W5). They could get into Crimea during the Bølling-Allerød interstadial along the Black Sea seashore from the Transcaucasia refugim, but the lineage vanished in the next Younger Dryas cooling. Younger samples (from Kr to BSN13 dated to c. 12000–2400 cal yr BP, in Fig. 2) belong to clade W7 and are grouped together with sequences from Serbia and Italy and pre- and post-LGM specimens from England. This suggests that Crimea was colonized from the Balkans in the beginning of Holocene. Noticeably, in the late Pleistocene and early Holocene, the Black Sea level was lower than today (Fig. 1). As a consequence, the Balkans and the Caucasus were connected to Crimea by land, which made the migration between these areas possible (Stanko, 2007; Demidenko, 2008). This is also supported by similarities in material culture from the Late Palaeolithic and Mesolithic sites in Bulgaria, Romania and Crimea, suggesting direct contacts among humans from these localities (Stanko, 2007). Despite more stable and favourable climatic conditions since the beginning of the Holocene, the post-LGM red deer population from clade W7 did not survive to modern times in Crimea. This population disappeared c. 2.5 ka and was not an ancestor of contemporary deer. Causes for this Holocene extinction are less clear. Based on soil development and pollen analysis, Cordova & Lehman (2005) reconstructed the environmental changes in Crimea during the Holocene. Their study revealed climate fluctuations and two events of dry and cool phases with a decrease in arboreal vegetation and spreading steppe taxa, dated to 11–7.5 and 4.2–3.5 ka, respectively. This climatic deterioration, however, is unlikely to have caused the extinction of the Holocene population, given that similar conditions are present today in some parts of Crimea, where red deer exist. Cordova & Lehman (2005) also found an increase in phryganoid plants, weeds and other cultivated species around 3.5 ka, suggesting a human impact on the Crimean environment. Between seventh and fifth centuries BC (2.65–2.45 ka), Greeks started to colonize Crimea. One of the earliest settlements on the northern Black Sea coast was established c. 700 yr BC (2.65 ka) on Berezan Island, which was then a part of the mainland (Dolukhanov & Shilik, 2007). The Greek colonization contributed not only to increasing amounts of pastoral plants but also to intensive deforestation, which took place c. 2.5 ka (Cordova & Lehman, 2005; Cordova, 2007). This human activity corresponds well with the observed Holocene red deer extinction, suggesting that intensive landscape management led to red deer habitat fragmentation and possibly to its extinction. Humans have been postulated as a cause of faunal extinction in many other studies (Burney et al., 2004; Burney & Flannery, 2005; Lorenzen et al., 2011; Rawlence et al., 2015). Some of these concentrated on the overkill hypothesis, while others suggested less direct reasons such as habitat alteration, biological invasion or diseases (Burney & Flannery, 2005). The second recolonization of the peninsula appears to have started ~1000 years after the post-LGM extinction of red deer, between the Late Antic and Mediaeval periods, as is supported by two Crimean specimens (SOKIL and PR2) dated to 1.2–0.9 cal ka (eighth to 11th centuries AD) – Figure 2. These, together with present samples, belong to Western haplogroup A (clade W9). This haplogroup originated on the Iberian Peninsula and spread across the whole of Europe (Skog et al., 2009; Niedziałkowska et al., 2011; Meiri et al., 2013). This suggests that Crimea was probably colonized naturally by immigrant red deer from Western Europe. However, translocations by humans cannot be excluded, as it is well known that this happened in Sardinia, Corsica, Crete and Ireland during the Neolithic period (Hajji et al., 2008; Carden et al., 2012; Harris, 2014). The reasons for such management are not well known, but studies by Sykes et al. (2006, 2011) and Davis & MacKinnon (2009) have showed that fallow deer (Dama dama) was introduced to Britain and Portugal by the Romans, probably for keeping them in special parks (vivaria), which enhanced their owner’s status. In the tenth century AD, Crimea was a part of the Eastern Roman Empire. Historical and archaeological studies suggest that deer was a totemic animal for Romans and a symbol of wealth and power. Moreover, hunting for red deer in Roman Britain was an expression of status (Allen et al., 2014). Therefore, a potential extinction of natural red deer population due to Greek activities in Crimea may have induced the need for its reintroduction. More studies on ancient red deer populations, especially from potential migration routes from Iberia to Crimea, are crucial to resolve the question about the origin of the contemporary red deer population in Crimea. Comparison of the 12 models with ABC analyses revealed that after the LGM, red deer likely colonized the Crimean Peninsula in two independent waves: at the end of the Pleistocene and in the Holocene. Those events were related to extinctions and replacements of populations. Similar events have been reported for the collared lemming, where population turnovers were probably caused by climate fluctuations and changes in vegetation (Brace et al., 2012; Palkopoulou et al., 2016). Red deer, however, has a wider environmental tolerance and is generally well adapted to various environments (Straus, 1981; Geist, 1998). Thus, climate as a single factor driving extensive changes in its population may be less likely than for small mammals. Correspondingly, red deer as an important game species has been hunted, translocated and introduced since ancient times (Hartl, Zachos & Nadlinger, 2003; Zachos & Hartl, 2011). We find no support for the hypothesis that Crimea was a Pleistocene refugium for red deer. At least two extinction events in the late Quaternary are in conflict with the classical refugium definition. Considering that our analyses were based on an mtDNA marker and thus reflect only the maternal lineage history, it would be interesting to also analyse nuclear markers to confirm our conclusions on red deer population history in Crimea. Nonetheless, the Crimean samples analysed here provide significant insights into the phylogeography and migration routes of red deer at a global scale. SUPPORTING INFORMATION Additional Supporting Information may be found in the online version of this article at the publisher’s web-site: Table S1. Detailed information about samples used in the study. Radiocarbon dates obtained in previous studies are marked by asterisks. GenBank accession numbers in italics indicate partial cytb sequences published by Stankovic et al. (2011). Table S2. Chronostratigraphy of layers in Emine-Bair-Khosar Cave. Table S3. Nucleotide substitution models for alignment sets used in three programs. In all models, we assumed five discrete categories in the case of gamma distributed rates. Table S4. Summary statistics used in ABC analysis. N – number of individuals; h – number of haplotypes; S – segregating sites; π – nucleotide diversity; PairDiff – pairwise differences between groups. Table S5. Model likelihoods and Bayes factors among the 12 models compared by approximate Bayesian computation. The Bayes factors have a row to column interpretation (the value at row i and column j is BFij = MLi/MLj). These values correspond to the analysis with a mutation rate of 10.0% per million years, after averaging the five smallest acceptance proportions (0.001 to 0.005% of the closest simulated datasets to the observed one). The best-fitted model was bolded. Table S6. Bayes factors among the 12 models compared by approximate Bayesian computation. The Bayes factors have a row to column interpretation (the value at row i and column j is BFij = MLi/MLj). These values correspond to the analysis with a mutation rate of 20.0% per million years, after averaging the ten smallest acceptance proportions (0.001–0.01% of the closest simulated datasets to the observed one). The number of acceptance proportions averaged for this result was higher than in tables with lower mutation rates due to many models presented a likelihood of 0.0 in the first five proportions. The best-fitted model was bolded. Table S7. Bayes factors among the 12 models compared by approximate Bayesian computation. The Bayes factors have a row to column interpretation (the value at row i and column j is BFij = MLi/MLj). These values correspond to the analysis with a mutation rate of 30.0% per million years, after averaging the ten smallest acceptance proportions (0.001 to 0.005% of the closest simulated datasets to the observed one). The number of acceptance proportions averaged for this result was higher than in tables with lower mutation rates due to several models presented a likelihood of 0.0 in the first five proportions. The best-fitted model was bolded. Table S8. Rate of error and statistical power in the model choice analysis by approximate Bayesian computation, as assessed by pseudo observed data sets (PODs). Figure S1. Cross section through layers in Emine-Bair-Khosar Cave. The position of red deer remains was marked. Uncalibrated radiocarbon dates of these samples and membership of DNA from the samples to phylogenetic clades were given. Figure S2. MrBayes tree based on the longer alignment of Cervini members with simplification of Eastern (clade E) and Western (clade W) lineages of red deer. The asterisk indicates that a given sequence is represented by more than one accession number. The shading of circles at nodes corresponds to the number of methods (out of five) which recovered the given node. Numbers at nodes, in the order shown, correspond to: posterior probabilities estimated in MrBayes (MB) and PhyloBayes (PB), bootstrap values obtained in PhyML(PH), TreeFinder (TF) and PAUP by neighbor-joining methods (NJ). Values of the posterior probabilities and bootstrap percentages ≤ 0.50 and 50%, respectively, were omitted or indicated by a dash ‘–’. Figure S3. Simplified maximum clade credibility tree obtained in BEAST for the longer alignment of Cervini members. Median ages expressed in ka BP are shown for all nodes. Sequences obtained from Crimea specimens are bolded. The Eastern (clade E) and Western (clade W) lineages of red deer were indicated. Figure S4. MrBayes tree based on the shorter alignment for the Eastern lineage of red deer. Clades E1-E3 represent C. canadensis (wapiti), whereas E4-E7 C. nippon (sika). Sequences newly obtained in this study from Crimea specimens are bolded. The asterisk indicates that a given sequence is represented by more than one accession number. The shading of circles at nodes corresponds to the number of methods (out of five) which recovered the given node. A probable age or calibrated dates are shown in square brackets. Numbers at nodes, in the order shown, correspond to: posterior probabilities estimated in MrBayes (MB) and PhyloBayes (PB), bootstrap values obtained in PhyML (PH), TreeFinder (TF) and PAUP by neighbor-joining methods (NJ). Values of the posterior probabilities and bootstrap percentages ≤ 0.50 and 50%, respectively, were omitted or indicated by a dash ‘–’. Figure S5. MrBayes tree based on the shorter alignment for the Western lineage of red deer. Sequences newly obtained in this study from Crimea specimens are bolded. The asterisk indicates that a given sequence is represented by more than one accession number. The shading of circles at nodes corresponds to the number of methods (out of five) which recovered the given node. A probable age or calibrated dates are shown in square brackets. Numbers at nodes, in the order shown, correspond to: posterior probabilities estimated in MrBayes (MB) and PhyloBayes (PB), bootstrap values obtained in PhyML (PH), TreeFinder (TF) and PAUP by neighbor joining methods (NJ). Values of the posterior probabilities and bootstrap percentages 0.50 and 50%, respectively, were omitted or indicated by a dash ‘–’. Box S1. Input file for scenario A with mutation rate 5.14% per million years for final simulations. Box S2. Input file for scenario B with mutation rate 5.14% per million years for final simulations. Box S3. Input file for scenario C with mutation rate 5.14% per million years for final simulations. Box S4. Input file for scenario D with mutation rate 5.14% per million years for final simulations. Box S5. Input file for scenario E with mutation rate 5.14% per million years for final simulations. Box S6. Input file for scenario F with mutation rate 5.14% per million years r for final simulations. Box S7. Input file for scenario G with mutation rate 5.14% per million years for final simulations. Box S8. Input file for scenario H with mutation rate 5.14% per million years for final simulations. Box S9. Input file for scenario I with mutation rate 5.14% per million years for final simulations. Box S10. Input file for scenario J with mutation rate 5.14% per million years for final simulations. Box S11. Input file for scenario K with mutation rate 5.14% per million years for final simulations. Box S12. Input file for scenario L with mutation rate 5.14% per million years for final simulations. ACKNOWLEDGEMENTS We are very grateful to three anonymous reviewers for their valuable comments and insightful remarks that have significantly improved the paper. We thank prof. Jan Burdukiewicz who supplied dates for prehistoric periods and archaeological cultures that were ascribed for some ancient samples. We are thankful to Dr R. Karasiewicz-Szypiorski for helpful discussions about Crimean history. The research for this paper was supported by Ministry of Science and Higher Education grant no. N N303 821140 and by the EU through the European Social Fund, contract number UDA-POKL.04.01.01-00-072/09-00. REFERENCES Allen MG , Baker K , Carden R , Madgwick R . 2014 . Chasing Sylvia’s Stag: placing deer in the countryside of Roman Britain . In: Baker K , Carden R , Madgwick R , eds. Deer and people . Oxbow, Oxford, UK: Windgather Press , 174 – 186 . Avdeienko Y . 2015 . The Last Interglacial vegetation and climate in the karstic areas of the Crimea and the Middle Dniester Basin . Physical Geography and Geomorphology (Kyiv) 4 : 102 – 109 . Banks W , Derrico F , Peterson A , Kageyama M , Colombeau G . 2008 . Reconstructing ecological niches and geographic distributions of caribou (Rangifer tarandus) and red deer (Cervus elaphus) during the Last Glacial Maximum . Quaternary Science Reviews 27 : 2568 – 2575 . Google Scholar CrossRef Search ADS Beaumont MA , Zhang W , Balding DJ . 2002 . Approximate Bayesian computation in population genetics . Genetics 162 : 2025 – 2035 . Bennett K , Provan J . 2008 . What do we mean by ‘refugia’ ? Quaternary Science Reviews 27 : 2449 – 2455 . Google Scholar CrossRef Search ADS Blackwell BA , Liang S , Golovanova LV , Doronichev VB , Skinner AR , Blickstein JI . 2005 . ESR at Treugol’naya Cave, Northern Caucasus Mt., Russia: dating Russia’s oldest archaeological site and paleoclimatic change in Oxygen Isotope Stage 11 . Applied Radiation and Isotopes 62 : 237 – 245 . Google Scholar CrossRef Search ADS Brace S , Palkopoulou E , Dalén L , Lister AM , Miller R , Otte M , Germonpré M , Blockley SPE , Stewart JR , Barnes I . 2012 . Serial population extinctions in a small mammal indicate Late Pleistocene ecosystem instability . Proceedings of the National Academy of Sciences of the United States of America 109 : 20532 – 20536 . Google Scholar CrossRef Search ADS Brugal JP , Croitor R . 2007 . Evolution, ecology and biochronology of herbivore associations in Europe during the last 3 million years . Quaternaire 18: 129 – 152 . Burney DA , Burney LP , Godfrey LR , Jungers WL , Goodman SM , Wright HT , Jull AJ . 2004 . A chronology for late prehistoric Madagascar . Journal of Human Evolution 47 : 25 – 63 . Google Scholar CrossRef Search ADS Burney DA , Flannery TF . 2005 . Fifty millennia of catastrophic extinctions after human contact . Trends in Ecology & Evolution 20 : 395 – 401 . Google Scholar CrossRef Search ADS Cameron RAD , Pokryszko BM , Horsák M . 2013 . Forest snail faunas from Crimea (Ukraine), an isolated and incomplete Pleistocene refugium . Biological Journal of the Linnean Society 109 : 424 – 433 . Google Scholar CrossRef Search ADS Carden RF , McDevitt AD , Zachos FE , Woodman PC , O’Toole P , Rose H , Monaghan NT , Campana MG , Bradley DG , Edwards CJ . 2012 . Phylogeographic, ancient DNA, fossil and morphometric analyses reveal ancient and modern introductions of a large mammal: the complex case of red deer (Cervus elaphus) in Ireland . Quaternary Science Reviews 42 : 74 – 84 . Google Scholar CrossRef Search ADS Clark PU , Dyke AS , Shakun JD , Carlson AE , Clark J , Wohlfarth B , Mitrovica JX , Hostetler SW , McCabe AM . 2009 . The Last Glacial Maximum . Science 325 : 710 – 714 . Google Scholar CrossRef Search ADS Cooper A , Turney C , Hughen KA , Barry W , Mcdonald HG , Bradshaw CJA . 2015 . Abrupt warming events drove Late Pleistocene Holarctic megafaunal turnover . Science Express 349 : 1 – 8 . Cordova CE . 2007 . Holocene Mediterranization of the southern Crimean vegetation: paleoecological records, regional climate change, and possible non-climatic influences . In: Yanko-Hombach V , Gilbert AS , Panin N , Dolukhanov PM , eds. The Black Sea flood question: changes in coastline, climate, and human settlement . Netherlands: Springer , 319 – 344 . Google Scholar CrossRef Search ADS Cordova CE , Lehman PH . 2005 . Holocene environmental change in southwestern Crimea (Ukraine) in pollen and soil records . The Holocene 15 : 263 – 277 . Google Scholar CrossRef Search ADS Criscuolo A . 2011 . morePhyML: improving the phylogenetic tree space exploration with PhyML 3 . Molecular Phylogenetics and Evolution 61 : 944 – 948 . Google Scholar CrossRef Search ADS Darriba D , Taboada GL , Doallo R , Posada D . 2012 . jModelTest 2: more models, new heuristics and parallel computing . Nature Methods 9 : 772 . Google Scholar CrossRef Search ADS Davis S , MacKinnon M . 2009 . Did the Romans bring fallow deer to Portugal ? Environmental Archaeology 14 : 15 – 26 . Google Scholar CrossRef Search ADS Demidenko YE . 2008 . The Early and Mid-Upper Palaeolithic of the North Black Sea region: an overview . Quartär 55 : 99 – 114 . Dolukhanov PM , Shilik KK . 2007 . Environment, sea-level changes, and human migrations in the northern Pontic area during late Pleistocene and Holocene times . In: Yanko-Hombach V , Gilbert AS , Panin N , Dolukhanov PM , eds. The Black Sea flood question: changes in coastline, climate, and human settlement . Netherlands: Springer , 297 – 318 . Google Scholar CrossRef Search ADS Drummond AJ , Suchard MA , Xie D , Rambaut A . 2012 . Bayesian phylogenetics with BEAUti and the BEAST 1.7 . Molecular Biology and Evolution 29 : 1969 – 1973 . Google Scholar CrossRef Search ADS Fagundes NJR , Ray N , Beaumont M , Neuenschwander S , Salzano FM , Bonatto SL , Excoffier L . 2007 . Statistical evaluation of alternative models of human evolution . Proceedings of the National Academy of Sciences of the United States of America 104 : 17614 – 17619 . Google Scholar CrossRef Search ADS Fernández-García JL , Carranza J , Martínez JG , Randi E . 2013 . Mitochondrial D-loop phylogeny signals two native Iberian red deer (Cervus elaphus) lineages genetically different to Western and Eastern European red deer and infers human-mediated translocations . Biodiversity and Conservation 23 : 537 – 554 . Google Scholar CrossRef Search ADS Gąsiorowski M , Hercman H , Ridush B , Stefaniak K . 2014 . Environment and climate of the Crimean Mountains during the Late Pleistocene inferred from stable isotope analysis of red deer (Cervus elaphus) bones from the Emine-Bair-Khosar Cave . Quaternary International 326–327: 243 – 249 . Geist V . 1998 . Deer of the world: their evolution, behaviour, and ecology . Mechanicsburg, PA: Stackpole Books . Gerasimenko N . 2007 . Environmental changes in the Crimean Mountains during the Last Interglacial–middle pleniglacial as recorded by pollen and lithopedology . Quaternary International 164–165 : 207 – 220 . Google Scholar CrossRef Search ADS Gilbert C , Ropiquet A , Hassanin A . 2006 . Mitochondrial and nuclear phylogenies of Cervidae (Mammalia, Ruminantia): systematics, morphology, and biogeography . Molecular Phylogenetics and Evolution 40 : 101 – 117 . Google Scholar CrossRef Search ADS Groves C . 2005 . The genus Cervus in eastern Eurasia . European Journal of Wildlife Research 52 : 14 – 22 . Google Scholar CrossRef Search ADS Guindon S , Dufayard JF , Lefort V , Anisimova M , Hordijk W , Gascuel O . 2010 . New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0 . Systematic Biology 59 : 307 – 321 . Google Scholar CrossRef Search ADS Hajji GM , Charfi-Cheikrouha F , Lorenzini R , Vigne JD , Hartl GB , Zachos FE . 2008 . Phylogeography and founder effect of the endangered Corsican red deer (Cervus elaphus corsicanus) . Biodiversity and Conservation 17 : 659 – 673 . Google Scholar CrossRef Search ADS Hall TA . 1999 . BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT . Nucleic Acids Symposium Series 41: 95 – 98 . Harris K . 2014 . Hunting, performance and incorporation: human–deer encounter in Late Bronze Age Crete . In: Baker K , Carden R , Madgwick R , eds. Deer and people . Oxbow, Oxford, UK: Windgather Press , 48 – 58 . Hartl GB , Zachos F , Nadlinger K . 2003 . Genetic diversity in European red deer (Cervus elaphus L.): anthropogenic influences on natural populations . Comptes Rendus Biologies 326 ( Suppl 1 ): S37 – S42 . Google Scholar CrossRef Search ADS Hassanin A , Delsuc F , Ropiquet A , Hammer C , Jansen van Vuuren B , Matthee C , Ruiz-Garcia M , Catzeflis F , Areskoug V , Nguyen TT , Couloux A . 2012 . Pattern and timing of diversification of Cetartiodactyla (Mammalia, Laurasiatheria), as revealed by a comprehensive analysis of mitochondrial genomes . Comptes Rendus Biologies 335 : 32 – 50 . Google Scholar CrossRef Search ADS Hewitt GM . 1996 . Some genetic consequences of ice ages, and their role, in divergence and speciation . Biological Journal of the Linnean Society 58 : 247 – 276 . Google Scholar CrossRef Search ADS Hewitt GM . 2004 . Genetic consequences of climatic oscillations in the Quaternary . Philosophical Transactions of the Royal Society of London B: Biological Sciences 359 : 183 – 195 ; discussion 195. Google Scholar CrossRef Search ADS Ho SY . 2014 . The changing face of the molecular evolutionary clock . Trends in Ecology & Evolution 29 : 496 – 503 . Google Scholar CrossRef Search ADS Hubberten H . 2004 . The periglacial climate and environment in northern Eurasia during the Last Glaciation . Quaternary Science Reviews 23 : 1333 – 1357 . Google Scholar CrossRef Search ADS Jobb G , von Haeseler A , Strimmer K . 2004 . TREEFINDER: a powerful graphical analysis environment for molecular phylogenetics . BMC Evolutionary Biology 4 : 18 . Google Scholar CrossRef Search ADS Krojerová-Prokešová J , Barančeková M , Koubek P . 2015 . Admixture of eastern and western European Red deer lineages as a result of postglacial recolonization of the Czech Republic (Central Europe) . The Journal of Heredity 106 : 375 – 385 . Google Scholar CrossRef Search ADS Kuwayama R , Ozawa T . 2000 . Phylogenetic relationships among European red deer, wapiti, and sika deer inferred from mitochondrial DNA sequences . Molecular Phylogenetics and Evolution 15 : 115 – 123 . Google Scholar CrossRef Search ADS Kuznetsova MV , Danilkin AA , Kholodova MV . 2012 . Phylogeography of red deer (Cervus elaphus): analysis of mtDNA cytochrome b polymorphism . Biology Bulletin 39 : 323 – 330 . Google Scholar CrossRef Search ADS Kuznetsova M , Volokh A , Domnich V , Tyshkevitch V , Danilkin A . 2007 . Molecular diversity of the red deer, Cervus elaphus (Cervidae), inhabiting East Europe . Vestnik Zoologii 41 : 505 – 509 . [in Russian] Lartillot N , Lepage T , Blanquart S . 2009 . PhyloBayes 3: a Bayesian software package for phylogenetic reconstruction and molecular dating . Bioinformatics 25 : 2286 – 2288 . Google Scholar CrossRef Search ADS Lorenzen ED , Nogués-Bravo D , Orlando L , Weinstock J , Binladen J , Marske KA , Ugan A , Borregaard MK , Gilbert MTP , Nielsen R , Ho SYW , Goebel T , Graf KE , Byers D , Stenderup JT , Rasmussen M , Campos PF , Leonard JA , Koepfli KP , Froese D , Zazula G , Stafford TW , Aaris-Sørensen K , Batra P , Haywood AM , Singarayer JS , Valdes PJ , Boeskorov G , Burns JA , Davydov SP , Haile J , Jenkins DL , Kosintsev P , Kuznetsova T , Lai X , Martin LD , McDonald HG , Mol D , Meldgaard M , Munch K , Stephan E , Sablin M , Sommer RS , Sipko T , Scott E , Suchard MA , Tikhonov A , Willerslev R , Wayne RK , Cooper A , Hofreiter M , Sher A , Shapiro B , Rahbek C , Willerslev E . 2011 . Species-specific responses of Late Quaternary megafauna to climate and humans . Nature 479 : 359 – 364 . Google Scholar CrossRef Search ADS Lorenzini R , Garofalo L . 2015 . Insights into the evolutionary history of Cervus (Cervidae, tribe Cervini) based on Bayesian analysis of mitochondrial marker sequences, with first indications for a new species . Journal of Zoological Systematics and Evolutionary Research 53 : 340 – 349 . Google Scholar CrossRef Search ADS Ludt CJ , Schroeder W , Rottmann O , Kuehn R . 2004 . Mitochondrial DNA phylogeography of red deer (Cervus elaphus) . Molecular Phylogenetics and Evolution 31 : 1064 – 1083 . Google Scholar CrossRef Search ADS Mahmut H , Masuda R , Onuma M , Takahashi M , Nagata J , Suzuki M , Ohtaishi N . 2002 . Molecular phylogeography of the red deer (Cervus elaphus) populations in Xinjiang of China: comparison with other Asian, European, and North American populations . Zoological Science 19 : 485 – 495 . Google Scholar CrossRef Search ADS Markov GG , Kuznetsova MV , Danilkin AA , Kholodova MV . 2012 . Analysis of genetic diversity of red deer (Cervus elaphus L.) in Bulgaria: implications for population conservation and sustainable management . Acta Zoologica Bulgarica 64 : 389 – 396 . Markov GG , Kuznetsova MV , Danilkin AA , Kholodova MV . 2015 . Genetic diversity of the red deer (Cervus elphus L.) in Hungary revealed by cytochrome b gene . Acta Zoologica Bulgarica 67 : 11 – 17 . Markova AK . 2011 . Small mammals from Palaeolithic sites of the Crimea . Quaternary International 231 : 22 – 27 . Google Scholar CrossRef Search ADS Markova AK , Simakova AN , Puzachenko AY . 2009 . Ecosystems of Eastern Europe at the time of maximum cooling of the Valdai glaciation (24–18 kyr BP) inferred from data on plant communities and mammal assemblages . Quaternary International 201 : 53 – 59 . Google Scholar CrossRef Search ADS Meiri M , Lister AM , Collins MJ , Tuross N , Goebel T , Blockley S , Zazula GD , van Doorn N , Dale Guthrie R , Boeskorov GG , Baryshnikov GF , Sher A , Barnes I . 2014 . Faunal record identifies Bering isthmus conditions as constraint to end-Pleistocene migration to the New World . Proceedings of the Royal Society of London B: Biological Sciences 281 : 20132167 . Google Scholar CrossRef Search ADS Meiri M , Lister AM , Higham TF , Stewart JR , Straus LG , Obermaier H , González Morales MR , Marín-Arroyo AB , Barnes I . 2013 . Late-glacial recolonization and phylogeography of European red deer (Cervus elaphus L.) . Molecular Ecology 22 : 4711 – 4722 . Google Scholar CrossRef Search ADS Nagata J , Masuda R , Tamate HB , Hamasaki Si , Ochiai K , Asada M , Tatsuzawa S , Suda K , Tado H , Yoshida MC . 1999 . Two genetically distinct lineages of the sika deer, Cervus nippon, in Japanese islands: comparison of mitochondrial D-loop region sequences . Molecular Phylogenetics and Evolution 13 : 511 – 519 . Google Scholar CrossRef Search ADS Niedziałkowska M , Jędrzejewska B , Honnen AC , Otto T , Sidorovich VE , Perzanowski K , Skog A , Hartl GB , Borowik T , Bunevich AN , Lang J , Zachos FE . 2011 . Molecular biogeography of red deer Cervus elaphus from eastern Europe: insights from mitochondrial DNA sequences . Acta Theriologica 56 : 1 – 12 . Google Scholar CrossRef Search ADS Palkopoulou E , Baca M , Abramson NI , Sablin M , Socha P , Nadachowski A , Prost S , Germonpré M , Kosintsev P , Smirnov NG , Vartanyan S , Ponomarev D , Nyström J , Nikolskiy P , Jass CN , Litvinov YN , Kalthoff DC , Grigoriev S , Fadeeva T , Douka A , Higham TF , Ersmark E , Pitulko V , Pavlova E , Stewart JR , Węgleński P , Stankovic A , Dalén L . 2016 . Synchronous genetic turnovers across Western Eurasia in Late Pleistocene collared lemmings . Global Change Biology 22 : 1710 – 1721 . doi:10.1111/gcb.13214 . Google Scholar CrossRef Search ADS Pérez-Espona S , Pérez-Barbería FJ , Goodall-Copestake WP , Jiggins CD , Gordon IJ , Pemberton JM . 2009 . Genetic diversity and population structure of Scottish Highland red deer (Cervus elaphus) populations: a mitochondrial survey . Heredity 102 : 199 – 210 . Google Scholar CrossRef Search ADS Pitra C , Fickel J , Meijaard E , Groves PC . 2004 . Evolution and phylogeny of old world deer . Molecular Phylogenetics and Evolution 33 : 880 – 895 . Google Scholar CrossRef Search ADS Polziehn RO , Strobeck C . 2002 . A phylogenetic comparison of red deer and wapiti using mitochondrial DNA . Molecular Phylogenetics and Evolution 22 : 342 – 356 . Google Scholar CrossRef Search ADS R Core Team . 2014 . R: a language and environment for statistical computing . Vienna : R Foundation for Statistical Computing . Rambaut A , Suchard M , Xie D , Drummond A . 2014 . Tracer v1. 6 . Ramsey CB , Scott EM , van der Plicht J . 2013 . Calibration for archaeological and environmental terrestrial samples in the time range 26–50 ka cal bp . Radiocarbon 55 : 2021 – 2027 . Google Scholar CrossRef Search ADS Randi E , Mucci N , Claro-Hergueta F , Bonnet A , Douzery EJP . 2001 . A mitochondrial DNA control region phylogeny of the Cervinae: speciation in Cervus and implications for conservation . Animal Conservation 4 : 1 – 11 . Google Scholar CrossRef Search ADS Rawlence NJ , Perry GLW , Smith IWG , Scofield RP , Tennyson AJD , Matisoo-Smith EA , Boessenkool S , Austin JJ , Waters JM . 2015 . Radiocarbon-dating and ancient DNA reveal rapid replacement of extinct prehistoric penguins . Quaternary Science Reviews 112 : 59 – 65 . Google Scholar CrossRef Search ADS Reimer PJ , Bard E , Bayliss A , Beck JW , Blackwell PG , Bronk Ramsey C , Buck CE , Cheng H , Edwards RL , Friedrich M . 2013 . IntCal13 and Marine13 radiocarbon age calibration curves 0–50,000 years cal BP . Radiocarbon 55 : 1869 – 1887 . Google Scholar CrossRef Search ADS Richter J . 2005 . Hasty foragers: the Crimea Island and Europe during the last interglacial . In: Chabai V , Richter J , Uthmeier T , eds. Kabazi II: Last Interglacial occupation, environment and subsistence, Palaeolithic sites of Crimea . Shlyakh: Simferopol-Cologne , 275 – 286 . Ridush B , Stefaniak K , Socha P , Proskurnyak Y , Marciszak A , Vremir M , Nadachowski A . 2013 . Emine-Bair-Khosar Cave in the Crimea, a huge bone accumulation of Late Pleistocene fauna . Quaternary International 284 : 151 – 160 . Google Scholar CrossRef Search ADS Ronquist F , Teslenko M , van der Mark P , Ayres DL , Darling A , Höhna S , Larget B , Liu L , Suchard MA , Huelsenbeck JP . 2012 . MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space . Systematic Biology 61 : 539 – 542 . Google Scholar CrossRef Search ADS Rosvold J , Røed KH , Hufthammer AK , Andersen R , Stenøien HK . 2012 . Reconstructing the history of a fragmented and heavily exploited red deer population using ancient and contemporary DNA . BMC Evolutionary Biology 12 : 191 . Google Scholar CrossRef Search ADS Sanderson MJ . 2002 . Estimating absolute rates of molecular evolution and divergence times: a penalized likelihood approach . Molecular Biology and Evolution 19 : 101 – 109 . Google Scholar CrossRef Search ADS Sandom C , Faurby S , Sandel B , Svenning JC . 2014 . Global late Quaternary megafauna extinctions linked to humans, not climate change . Proceedings of the Royal Society of London B: Biological Sciences 281 : 20133254 . Google Scholar CrossRef Search ADS Sandoval-Castellanos E , Palkopoulou E , Dalén L . 2014 . Back to BaySICS: a user-friendly program for Bayesian Statistical Inference from Coalescent Simulations . PLoS ONE 9 : e98011 . Google Scholar CrossRef Search ADS Skog A , Zachos FE , Rueness EK , Feulner PGD , Mysterud A , Langvatn R , Lorenzini R , Hmwe SS , Lehoczky I , Hartl GB , Stenseth NC , Jakobsen KS . 2009 . Phylogeography of red deer (Cervus elaphus) in Europe . Journal of Biogeography 36 : 66 – 77 . Google Scholar CrossRef Search ADS Sommer RS , Zachos FE . 2009 . Fossil evidence and phylogeography of temperate species: ‘glacial refugia’ and post-glacial recolonization . Journal of Biogeography 36 : 2013 – 2020 . Google Scholar CrossRef Search ADS Sommer RS , Zachos FE , Street M , Jöris O , Skog A , Benecke N . 2008 . Late Quaternary distribution dynamics and phylogeography of the red deer (Cervus elaphus) in Europe . Quaternary Science Reviews 27 : 714 – 733 . Google Scholar CrossRef Search ADS Stanko VN . 2007 . Fluctuations in the level of the Black Sea and Mesolithic settlement of the northern Pontic area . In: Yanko-Hombach V , Gilbert AS , Panin N , Dolukhanov PM , eds. The Black Sea flood question: changes in coastline, climate, and human settlement . Netherlands: Springer , 371 – 385 . Google Scholar CrossRef Search ADS Stankovic A , Doan K , Mackiewicz P , Ridush B , Baca M , Gromadka R , Socha P , Weglenski P , Nadachowski A , Stefaniak K . 2011 . First ancient DNA sequences of the Late Pleistocene red deer (Cervus elaphus) from the Crimea, Ukraine . Quaternary International 245 : 262 – 267 . Google Scholar CrossRef Search ADS Stewart JR , Lister AM , Barnes I , Dalén L . 2010 . Refugia revisited: individualistic responses of species in space and time . Proceedings of the Royal Society of London B: Biological Sciences 277 : 661 – 71 . Google Scholar CrossRef Search ADS Stiller M , Knapp M , Stenzel U , Hofreiter M , Meyer M . 2009 . Direct multiplex sequencing (DMPS) – a novel method for targeted high-throughput sequencing of ancient and highly degraded DNA . Genome Research 19 : 1843 – 1848 . Google Scholar CrossRef Search ADS Straus LG . 1981 . On the habitat and diet of Cervus elaphus . Munibe 33 : 175 – 182 . Swofford DL . 1998 . Phylogenetic analysis using parsimony (* and other methods), Version 4 . Sykes NJ , Baker KH , Carden RF , Higham TFG , Hoelzel AR , Stevens RE . 2011 . New evidence for the establishment and management of the European fallow deer (Dama dama dama) in Roman Britain . Journal of Archaeological Science 38 : 156 – 165 . Google Scholar CrossRef Search ADS Sykes N , White J , Hayes T , Palmer M . 2006 . Tracking animals using strontium isotopes in teeth: the role of fallow deer (Dama dama) in Roman Britain . 80 : 948 – 960 . Taberlet P , Fumagalli L , Wust-Saucy AG , Cosson JF . 1998 . Comparative phylogeography and postglacial colonization routes in Europe . Molecular Ecology 7 : 453 – 464 . Google Scholar CrossRef Search ADS Tamura K , Stecher G , Peterson D , Filipski A , Kumar S . 2013 . MEGA6: Molecular Evolutionary Genetics Analysis version 6.0 . Molecular Biology and Evolution 30 : 2725 – 2729 . Google Scholar CrossRef Search ADS Vremir MM , Ridush B . 2005 . The Emine-Bair-Khosar ‘Mega Trap’ (Ukraine) . Mitteilungen der Kommission für Quartärforschung Österreichischen Akademie der Wissenschaften 14: 235 – 239 . Zachos FE , Hartl GB . 2011 . Phylogeography, population genetics and conservation of the European red deer Cervus elaphus . Mammal Review 41 : 138 – 150 . Google Scholar CrossRef Search ADS © 2017 The Linnean Society of London, Zoological Journal of the Linnean Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Zoological Journal of the Linnean Society Oxford University Press

The history of Crimean red deer population and Cervus phylogeography in Eurasia

Loading next page...
 
/lp/ou_press/the-history-of-crimean-red-deer-population-and-cervus-phylogeography-MeEt3uQpAW
Publisher
Oxford University Press
Copyright
© 2017 The Linnean Society of London, Zoological Journal of the Linnean Society
ISSN
0024-4082
eISSN
1096-3642
D.O.I.
10.1093/zoolinnean/zlx065
Publisher site
See Article on Publisher Site

Abstract

Abstract The present distribution of many species is a result of climatic changes during the Pleistocene and human activity. The impact of climate has been accompanied by restrictions of populations into refugia during glacial periods, and subsequent expansions during more favourable conditions, whereas human influence has been associated with hunting practices and translocations. One mammalian species that has been subject to such transformations is the red deer, Cervus elaphus, but the exact nature of these changes has been difficult to determine using only modern DNA. In this study, we obtained new cytochrome b sequences from subfossil remains of deer found in the Crimean Peninsula. A comparison of these sequences with the available recent and ancient sequences allowed to us to reconstruct phylogeographic relationships between Cervus lineages and to determine their potential migration routes at both local and Eurasian scales. Our analyses showed that the Crimean Peninsula was not a glacial refugium for red deer, but rather that red deer colonized Crimea in three independent waves from both Western and Eastern red deer populations. The immigrations were related to local extinctions and replacements of native populations. ancient DNA, Cervus elaphus, Crimea, phylogeny, phylogeography, Pleistocene, red deer INTRODUCTION During the Late Pleistocene, both environmental change and human activity influenced the distribution and genetic diversity of many species (Hewitt, 1996; Burney & Flannery, 2005; Stewart et al., 2010; Lorenzen et al., 2011; Sandom et al., 2014; Cooper et al., 2015). These factors led to the extinction of many taxa and shaped the modern-day phylogeographic patterns of those that survived. Genetic analyses suggested repeated retreats of ancient populations into refugia during unfavourable conditions and subsequent expansions after a suitable environment stabilized (Hewitt, 1996; Taberlet et al., 1998; Bennett & Provan, 2008). One of the species subjected to such climatic and environmental change was the red deer (Cervus elaphus). Its current distribution is assumed to be strongly influenced by its colonization history during the Late Pleistocene and early Holocene and by human activities throughout the Holocene. Phylogeographic analyses showed two main clades of this species (Kuwayama & Ozawa, 2000; Ludt et al., 2004; Skog et al., 2009). Individuals in Asia and North America have been assigned to the Eastern clade, and European, Middle Eastern and African populations belong to the Western clade. These clades represent populations that are adapted to different climatic conditions (Geist, 1998). The European red deer has a mixed feeding strategy, whereas the Eastern deer is a cold-adapted open-country grazer (Geist, 1998). Estimates of divergence times among the Western and Eastern lineages differ among studies and range from less than a million (Kuwayama & Ozawa, 2000; Mahmut et al., 2002; Polziehn & Strobeck, 2002) to 6–7 Mya (Randi et al., 2001; Ludt et al., 2004). Phylogeographic analyses of European red deer populations have shown that their isolation in Pleistocene refugia had a major impact on their genetic diversity and distribution (Skog et al., 2009; Niedziałkowska et al., 2011). These studies found three mtDNA lineages, marked as haplogroups A, B and C. It has been inferred that these haplogroups correspond to isolation in Iberian, African/Sardinian and Balkan refugia, respectively. Nowadays, haplogroup A is present in Western, Northern and Central Europe. Haplogroup B is restricted to North Africa and Sardinia, whereas haplogroup C is found in Southern and Eastern Europe. Human management over the last centuries is believed to have disrupted the genetic diversity of red deer populations in Europe that arose during the long-term isolation and subsequent post-glacial decolonization. Niedziałkowska et al. (2011) found that haplogroups from lineages A and C occur together in Poland, Belarus and Lithuania. This could be a result of either natural overlap and contact zones or human-mediated translocations that occurred at the end of the 19th and beginning of the 20th centuries. A similar mixture of mtDNA lineages has been found in populations from Bulgaria, Hungary and the Czech Republic (Markov et al., 2012, 2015; Krojerová-Prokešová, Barančeková & Koubek, 2015). Red deer translocations have also been suggested for populations in Germany (Ludt et al., 2004), Scotland (Carden et al., 2012) and Spain (Fernández-García et al., 2013). Moreover, Kuznetsova et al. (2007) found the mtDNA lineages of both Western and Eastern red deer in the Ukrainian population, which probably resulted from the introduction of Cervus elaphus maral to Europe. Therefore, to better understand the impact of climate and environmental conditions on the natural distribution of red deer populations during in the last glaciation, it is crucial to analyse samples from the period before the intensive human management. This can be achieved by ancient DNA (aDNA) analyses, as demonstrated in earlier studies of red deer remains, which indicated that the Late Quaternary history of red deer was more complex than previously assumed (Stankovic et al., 2011; Meiri et al., 2013). Its past genetic diversity was greater than at present, and the distribution of mtDNA lineages was also different. Meiri et al. (2013) noticed that mtDNA lineages only found today in Southern and Eastern Europe were also present in England and Spain before the Last Glacial Maximum (LGM). Analyses of Crimean red deer (Stankovic et al., 2011) indicated that both Western and Eastern lineages were present on this peninsula in the Late Pleistocene, suggesting that conclusions about isolation by distance, which were based on contemporary samples, could be misleading. The Crimean Peninsula has been proposed as a Pleistocene refugium for several species. Biome reconstructions have revealed that forest and steppe-forest vegetation with pine and broadleaved species persisted in Crimea through the last glaciation (Markova, Simakova & Puzachenko, 2009), although Cameron, Pokryszko & Horsák (2013) have argued that full forest cover probably did not survive. Environmental conditions appear to have been stable enough to create a refugium for small mammals, snails and trees (Gerasimenko, 2007; Markova et al., 2009; Markova, 2011; Cameron et al., 2013). Black Sea-level fluctuations in the Pleistocene may periodically have cut off the peninsula from the mainland. The sea level during the Eemian interglacial was few metres higher than today, which could result in an island (Richter, 2005). Crimea became the part of the mainland during colder periods due to sea-level decrease. This level was lower by 30–120 m than present during the Vistulian Interpleniglacial (MIS3), when Crimea was connected with Dobruja, a historical area located in the eastern part of Romania and Bulgaria (Stanko, 2007; Demidenko, 2008). The sea-level fluctuations may therefore have influenced the history of the red deer in the region. Deer populations may have been isolated in the peninsula during warmer events and migrations from adjacent regions may have occurred during colder phases. Similar landscape changes did not happen in other European refugia, which makes Crimea an unusual and interesting area for studying species’ response to Pleistocene climate and geographic changes. In order to test the hypothesis by Stankovic et al. (2011) that the Crimean Peninsula served as Europe’s refugium, from which red deer colonized the Eurasian continent, we obtained additional aDNA from red deer specimens dated to the Late Pleistocene and Holocene. If this hypothesis is correct, we would expect that the Late Pleistocene red deer population would survive the LGM. Therefore, the genetic continuity between samples from before and after the LGM should be observed. Since ancient samples from Crimea represented both Western and Eastern lineages of red deer (Stankovic et al., 2011), we also included all Eurasian available red deer sequences published so far, to reconstruct phylogeographic relationships and determine their potential migration routes at a Eurasian scale. MATERIAL AND METHODS Samples and DNA analyses We analysed 33 samples from red deer remains excavated from ten caves located on three plateaus in the Crimean Mountains (Fig. 1; Supporting Information, Table S1). The majority of the samples were obtained from the Emine-Bair-Khosar Cave (EBK), which is an important palaeontological site rich in vertebrate remains (Vremir & Ridush, 2005; Ridush et al., 2013) (Supporting Information, Fig. S1 and Table S2). Seventeen samples yielding mtDNA were radiocarbon dated in the Poznań Radiocarbon Laboratory and Gliwice Absolute Dating Methods Centre in Poland. Five other samples had been dated previously. The age of the samples ranged from c. 47000 to 900 cal yr BP (Fig. 2; Supporting Information, Table S1), and thus covered a broad time range from the Late Pleistocene to the Holocene. Figure 1. View largeDownload slide A, Crimean Peninsula with the research area (in rectangle) and palaeocoastlines in four time periods. B, research area with the location of caves: 1 – Emine-Bair-Khosar, 2 – K18, 3 – Pavucha, 4 – Primula, 5 – Binbash Koba, 6 – Krapivnyi Grotto, 7 – Dizel, 8 – Bilosnizhka. 9 – Lavrushka, 10 – Sokil. Figure 1. View largeDownload slide A, Crimean Peninsula with the research area (in rectangle) and palaeocoastlines in four time periods. B, research area with the location of caves: 1 – Emine-Bair-Khosar, 2 – K18, 3 – Pavucha, 4 – Primula, 5 – Binbash Koba, 6 – Krapivnyi Grotto, 7 – Dizel, 8 – Bilosnizhka. 9 – Lavrushka, 10 – Sokil. Figure 2. View largeDownload slide Distribution of calibrated radiocarbon dates in time for red deer samples from Crimea. Median (circle) and 95.4% confidence interval (whiskers) are shown. The circles with arrows mean open dates. Phylogenetic clades, to which the samples belong, are shown in parentheses. Figure 2. View largeDownload slide Distribution of calibrated radiocarbon dates in time for red deer samples from Crimea. Median (circle) and 95.4% confidence interval (whiskers) are shown. The circles with arrows mean open dates. Phylogenetic clades, to which the samples belong, are shown in parentheses. The DNA extraction was performed by the Institute of Genetics and Biotechnology of the University of Warsaw, in a laboratory dedicated to aDNA analyses. The amplification of the cytochrome b (cytb) sequences (1140 bp) was performed in multiplex PCR using 12 overlapping primer pairs described in Stankovic et al. (2011). Full details regarding DNA extraction and amplification are given in Supporting Information. Multiplex PCR products were used to prepare libraries for pyrosequencing following the protocol from Stiller et al. (2009) and were sequenced on a Roche 454 sequencing platform at the Institute of Biochemistry and Biophysics, PAS. The resulting sequences were trimmed in Bioedit v. 7.2.5 (Hall, 1999) and assembled into contigs using SeqMan Pro software (DNASTAR, Inc.). To ensure authenticity of the sequence data, each sample was amplified in two independent multiplex PCRs and contigs from both replicates were used to create consensus sequences according to guidelines proposed by Stiller et al. (2009). Phylogenetic analyses Data sets were obtained by comprehensive searches of the GenBank database for cytochrome b homologs in members of the Cervini tribe. Many of the red deer sequences studied here, including fossil samples, were obtained by various research groups, such as Ludt et al. (2004), Skog et al. (2009), Kuznetsova, Danilkin & Kholodova (2012) and Meiri et al. (2013, 2014). Phylogenetic analyses were based on two alignments comprising non-redundant sequences: a longer set of 190 sequences, 1140 bp in length and a shorter set of 175 sequences with 423 bp. We applied five approaches for phylogenetic reconstruction: Bayesian inference in MrBayes (Ronquist et al., 2012) and PhyloBayes (Lartillot, Lepage & Blanquart, 2009), maximum likelihood analyses in TreeFinder (Jobb, von Haeseler & Strimmer, 2004) and morePhyML (Criscuolo, 2011) based on PhyML (Guindon et al., 2010), and the neighbor-joining method in PAUP* (Swofford, 1998). In the MrBayes analyses, we assumed separate mixed + I + Γ(5) models for each of three codon positions and two independent runs, each using 32 and 16 Markov chains for the longer and shorter sets, respectively. Trees were sampled every 100 generations for 10000000 generations. In the final analysis, we selected trees from the last 3062000 (the longer set) and 4664000 (the shorter set) generations that reached convergence. In the PhyloBayes approach, two independent Markov chains were run through 100000 (the longer set) and 200000 (the shorter set) cycles assuming the CAT GTR + Γ(5) model. After obtaining convergence, the last 50000 and 10000 trees from each chain were collected to compute the posterior consensus, respectively. In TreeFinder, we applied a search depth set to 2, separate substitution models for each of three codon positions according to the TreeFinder Propose Model (Supporting Information, Table S3). Trees inferred with morePhyML and PAUP were based on the best-fit substitution models found in jModeltest (Darriba et al., 2012) among 1624 candidate models (Supporting Information, Table S3). The best search algorithm NNI + SPR was applied in morePhyML. To assess the significance of branches, non-parametric bootstrap analyses were performed on 1000 replicates in TreeFinder, PhyML and PAUP. Molecular dating Divergence times were estimated with the software BEAST (Drummond et al., 2012) on the full-length alignment comprising 216 sequences, using as many samples as possible with assigned dates and geographic localizations. We used a fixed MrBayes tree that was inferred according to the procedure described above, but with 32 Markov chains and the last 3861000 generations to summarize the final tree. Based on that, we created an ultrametric tree using Sanderson’s approach (Sanderson, 2002) in the chronopl function in R (R Core Team, 2014). We applied separate substitution models for each of three codon positions selected in jModeltest among 1624 candidate models (Supporting Information, Table S2). The age of the root of the tree was set as a uniform prior with the lower and upper limits 5.1 and 7.3 Mya. This age corresponds to the separation time of Axis and Rucervus from other Cervini based on three methods (Hassanin et al., 2012). Moreover, we used 25 tip dates associated with sequences obtained from fossil and historical samples. To unify the dates, we applied the same calibration procedure in OxCal 4.2.3 (Ramsey, Scott & van der Plicht, 2013) using curve IntCal 13 (Reimer et al., 2013) for original uncalibrated dates. We tested strict and uncorrelated lognormal relaxed clock models with a coalescent constant size tree and separate clock models for particular data partitions. Finally, we selected the relaxed clock model because the SD of the relaxed clock revealed significant variation among branches. Posterior distributions of parameters were estimated for 50000000 generations with a sampling frequency of 1000 steps. The convergence and sufficient sampling were checked using Tracer (Rambaut et al., 2014). The phylogenetic trees were summarized in TreeAnnotator (Drummond et al., 2012) with a 10% burn-in. Approximate Bayesian computation Based on archaeological and historical data, we defined a set of the 12 most likely models that described the temporal pattern in genetic diversity of Crimean red deer (Fig. 3). The best-fitted model was selected by approximate Bayesian computation (ABC) (Beaumont, Zhang & Balding, 2002). For defining the models, we considered the following scenarios: Figure 3. View largeDownload slide Population models considered in approximate Bayesian computation analyses. The most likely model H assuming 5.14 and 10% mutation rate per Myr was outlined by the red rectangle. Figure 3. View largeDownload slide Population models considered in approximate Bayesian computation analyses. The most likely model H assuming 5.14 and 10% mutation rate per Myr was outlined by the red rectangle. Two possible scenarios regarding the pre-LGM origin of Crimean population: one including a long-standing panmictic population (Fig. 3A, C, E, G, I, K) and the second assuming an admixture of Western and Eastern populations (Fig. 3B, D, F, H, J, L); Two possible scenarios regarding the post-LGM (until c. 2000 cal yr BP) origin of the Crimean population: in one scenario, Crimea was populated by deer that survived locally (Fig. 3A–F). In the other, Crimea was colonized by Western European lineages (Fig. 3H–L); Three scenarios regarding the present origin of the Crimean population: in one scenario, the present population was derived from the population inhabiting Crimea before 2000 cal yr BP (Fig. 3A, B, K, L); in the second scenario, the present population split from a common ancestor with the population that colonized Crimea after the LGM and until 2000 cal yr BP. Here we considered two split times: before the LGM (Fig. 3I, J) and after the LGM (Fig. 3C, D). In the third scenario, the present population was colonized from Western European lineages that were unrelated to the post-LGM (till c. 2000 cal yr BP) population (Fig. 3E–H). The 12 tested models resulted from combining these scenarios. Simulations were performed under a coalescent framework, where unknown parameters and their uncertainness were assessed by Bayesian prior probability distributions. We carried out two rounds of pilot simulations to optimize the parameters. The priors for effective population size (Ne) had an exponential distribution in pilot simulations to better screen an ample range. Then, we replaced them with uniform distributions in the final simulations. The times to population splits or demographic changes were set by uniform priors as well. Uncertainty in radiocarbon dating of sample ages was included by assigning normal priors (except in four cases where exponential priors were used because of infinite radiocarbon dates; see Supporting Information). We employed a 4-year generation time, based on the literature (Pérez-Espona et al., 2009; Rosvold et al., 2012) and a transition/transversion (Tr/Ts) bias of 0.945, a gamma shape parameter of 0.05 and a mutation rate of 5.14% per Myr (Skog et al., 2009). We calculated the basic parameters: the number of haplotypes, private haplotypes and segregating sites and nucleotide diversity for three sets of samples grouped by time periods: pre-LGM, post-LGM to 2000 cal yr BP and the present. In addition, we employed the number of pairwise differences and FST among pairs of consecutive samples (Supporting Information, Table S4). We ran 2000000 simulations for each model in the final analyses. The best-fitting model was selected with Bayes factors (BFs) that were calculated as ratios of model likelihoods between pairs of models for all 66 pairwise comparisons of the 12 models (Supporting Information, Table S5). The model likelihoods were estimated by (1) a direct approach where they were directly equalled to their relative acceptance ratios, (2) a weighted approach, which also employs the acceptance ratios after being weighted with an Epanechnikov kernel evaluated in the Euclidean distances to the observed data and (3) a logistic regression in which summary statistics represent explanatory variables of model likelihoods (Fagundes et al., 2007). The consistency of the model likelihoods and the BFs was assessed by applying the procedure with 10 different acceptance proportions (from 0.001 to 0.01%). Since the mutation rate, and possible dependency on temporal scale (Ho, 2014), could impact our estimated parameters and model choice, we tested the robustness of our results with a range of different mutation rates: 10, 20 and 30% per Myr. Finally, to assess the probability of correct model selection, we estimated the statistical power of the analyses with pseudo observed data sets (PODs), where simulated data sets play the role of real ones. We generated 1000 PODs, which were analysed with the same reference tables and 2000000 simulations each. Simulations, model choice and POD analyses were performed in the software BaySICS (Sandoval-Castellanos, Palkopoulou & Dalén, 2014), whereas the Tr/Ts bias and gamma parameter were calculated in MEGA 6 (Tamura et al., 2013). RESULTS Phylogenetic relationships within red deer lineages including Crimean deer sequences We successfully amplified and obtained cytochrome b sequence from 24 samples, which gave us a 73% success rate of DNA extraction for the subfossil samples. We identified 16 haplotypes. None of the negative controls yielded DNA. The sequences were aligned to all available cytochrome b sequences from Cervini members including ancient and modern samples of C. elaphus. The five applied approaches produced similar trees. Deer sequences assigned to the Cervus genus do not form a monophyletic clade but are split into Western and Eastern subclades, with strong and moderate support, respectively (Supporting Information, Figs S2, S3). The Eastern subclade is very significantly grouped with the Thorold’s deer Przewalskium albirostris sequences. All newly obtained sequences from the Crimean samples group significantly either with the Western or Eastern lineages (detailed phylogenies in Figs 4, 5; Supporting Information, Figs S4, S5). To study the detailed relationships between Cervus sequences, we considered the shorter set (Supporting Information, Figs S4, S5) in addition to the longer one (Figs 4, 5) because it contained important subfossil samples. Figure 4. View largeDownload slide MrBayes tree based on the longer alignment for the Eastern lineage of red deer (the clade E in Supporting information, Fig. S2). Clades E2 and E3 represent C. elaphus canadensis (wapiti), whereas E4–E7 C. elaphus nippon (sika). Sequences newly obtained in this study from Crimea specimens are in bold. The asterisk indicates that a given sequence is represented by more than one accession number. The shaded circles at nodes correspond to the number of methods (out of five) that recovered the given node. Probable age, time period or calibrated dates of samples are shown in square brackets. Numbers at nodes, in the order shown, correspond to posterior probabilities estimated in MrBayes (MB) and PhyloBayes (PB), bootstrap values obtained in PhyML (PH), TreeFinder (TF) and PAUP by the neighbor-joining (NJ) method. Values of the posterior probabilities and bootstrap percentages ≤ 0.50 and 50%, respectively, were omitted or indicated by a dash ‘–’. Figure 4. View largeDownload slide MrBayes tree based on the longer alignment for the Eastern lineage of red deer (the clade E in Supporting information, Fig. S2). Clades E2 and E3 represent C. elaphus canadensis (wapiti), whereas E4–E7 C. elaphus nippon (sika). Sequences newly obtained in this study from Crimea specimens are in bold. The asterisk indicates that a given sequence is represented by more than one accession number. The shaded circles at nodes correspond to the number of methods (out of five) that recovered the given node. Probable age, time period or calibrated dates of samples are shown in square brackets. Numbers at nodes, in the order shown, correspond to posterior probabilities estimated in MrBayes (MB) and PhyloBayes (PB), bootstrap values obtained in PhyML (PH), TreeFinder (TF) and PAUP by the neighbor-joining (NJ) method. Values of the posterior probabilities and bootstrap percentages ≤ 0.50 and 50%, respectively, were omitted or indicated by a dash ‘–’. Figure 5. View largeDownload slide MrBayes tree based on the longer alignment for the Western lineage of red deer (the clade W in Supporting Information, Fig. S2). Sequences newly obtained in this study from Crimean specimens are shown in bold. For further explanation, see Figure 4. Figure 5. View largeDownload slide MrBayes tree based on the longer alignment for the Western lineage of red deer (the clade W in Supporting Information, Fig. S2). Sequences newly obtained in this study from Crimean specimens are shown in bold. For further explanation, see Figure 4. In the Eastern lineage, two main well-supported groups representing wapiti, Cervus canadensis, and sika deer, Cervus nippon, can be recognized (Fig. 4; Supporting Information, Fig. S4). The wapiti group consists of clades marked as E1 (only in the tree for the smaller set, Supporting Information, Fig. S4), E2 and E3. The basal clade E1 comprises exclusively ancient sequences from Alaska and north-eastern Siberia (Meiri et al., 2014). The pre-LGM samples from north-eastern Siberia (Meiri et al., 2014) are located in the E2 and E3 clades too. The E2 clade also includes other samples from Central and central-eastern Asia. The deer samples in the E3 clade show a much wider distribution. Besides Central Asia, north-eastern Siberia and South Korea, they were also found in North America and Crimea, where they are represented by only pre-LGM samples. In the tree for the longer set, all Crimean sequences branch off at the base of E3 clade, but for the shorter set, they are placed inside the clade without any support (Fig. 4; Supporting Information, Fig. S4). In the sika deer group, four clades can be recognized (Fig. 4; Supporting Information, Fig. S4) including Asian (clades E4 and E5) and Japanese deer separated into its ‘northern’ (clade E6) and ‘southern’ (clade E7) (Nagata et al., 1999; Groves, 2005). Deer introduced to the Czech Republic are also present in these clades. However, no Crimean samples were reported among the sika deer clade. The Western red deer lineage is diversified into nine distinct clades, which are inferred by at least three out of the five methods applied (Fig. 5; Supporting Information, Fig. S5). Sequences from Crimean deer are distributed among four clades (W5, W6, W7 and W9). The Late Pleistocene samples from Crimea (post-LGM) are present in clade W5, which also contains other sequences from the southern part of Eastern Europe: modern samples from Krasnodarsky Krai, Voronezh and the Caucasus (from the pre-LGM period) (Meiri et al., 2013). Much older Crimea sequences from the pre-LGM are grouped with sequences from the Balkans, Italy, Central Europe and Turkey in the W6 clade (corresponding to haplogroup C). In the smaller set tree (Supporting Information, Fig. S5), this clade also includes pre-LGM specimens from the Caucasus and the Urals (Meiri et al., 2013) as well as a modern Crimean sequence. In the longer alignment tree (Fig. 5), the W7 clade contains exclusively post-LGM Crimean sequences, whereas in the shorter set tree (Supporting Information, Fig. S5), it includes additionally a modern sequence from Italy, two post-LGM sequences from Serbia and two ancient sequences from England, with pre- and post-LGM ages (Meiri et al., 2013), respectively. The W9 clade (corresponding to haplogroup A) has the highest sequence diversity. Besides Mediaeval (eighth to 11th centuries) and modern sequences from Crimea, it comprises deer mainly from Western, Central and Northern Europe including Great Britain, and sequences from regions far to the east, in Russian Karelia and Belarus, and deer introduced to South Korea and New Zealand. In the tree based on the short alignment (Supporting Information, Fig. S5), the basal position of this clade is occupied by a pre-LGM sample from Belgium (Meiri et al., 2013). Post-LGM specimens from Spain are branched off next. The deer from Crimea was not found in the basal clades of the Western lineage: W1 including Cervus elaphus yarkandensis and bactrianus from China and Uzbekistan; W2 (in the shorter alignment tree only, Supporting Information, Fig. S5) containing exclusively post-LGM sequences from Spain (Meiri et al., 2013); W3 (in the shorter alignment tree only, Supporting Information, Fig. S5) with a historical Jordan sample (Meiri et al., 2013); and W4 comprising the sequences of C. elaphus maral from Turkey and Iran. The W8 clade (corresponding to haplogroup B) is also absent from Crimean samples and includes deer from Tunisia, Sardinia, Hungary, the Swiss Alps, a deer introduced to Germany and pre-LGM samples from the Caucasus – Supporting Information, Figure S5. Population history of Crimean deer In the ABC model choice analysis, the same model H was selected for mutation rates of 5.14 and 10% per Myr (Fig. 3H; Table 1; Supporting Information, Table S5). This model assumed that the pre-LGM population was an admixture of Eastern (Asian) and Western (Europe) lineages, whereas the two post-LGM populations (post-LGM to 2000 yr BP, and present) resulted from two unrelated waves of colonization that originated in the West before 50 ka. The BF of the chosen model was relatively large (Table 1; Supporting Information, Table S5). For mutation rates of 20 and 30% per Myr, model G instead obtained the highest support (Supporting Information, Tables S6 and S7). The model G assumed that the pre-LGM Crimean population was panmictic and long standing. Table 1. Model likelihoods and Bayes factors among the 12 models compared by approximate Bayesian computation Model Model likelihoods A B C D E F G H I J K L A 0.0006 – 0.348 0.071 0.128 0.021 0.021 0.004 0.001 0.006 0.003 0.095 0.034 B 0.0018 2.871 – 0.205 0.368 0.059 0.059 0.012 0.004 0.017 0.009 0.272 0.097 C 0.0087 14.000 4.876 – 1.793 0.289 0.288 0.056 0.020 0.085 0.044 1.327 0.474 D 0.0049 7.806 2.719 0.558 – 0.161 0.161 0.031 0.011 0.048 0.024 0.740 0.264 E 0.0300 48.419 16.865 3.459 6.202 – 0.996 0.194 0.068 0.295 0.151 4.590 1.640 F 0.0302 48.613 16.933 3.472 6.227 1.004 – 0.195 0.068 0.296 0.152 4.609 1.647 G 0.1547 249.355 86.854 17.811 31.942 5.150 5.129 – 0.348 1.519 0.779 23.639 8.448 H 0.4439 715.903 249.360 51.136 91.707 14.785 14.727 2.871 – 4.362 2.237 67.869 24.255 I 0.1018 164.129 57.169 11.724 21.025 3.390 3.376 0.658 0.229 – 0.513 15.560 5.561 J 0.1985 320.032 111.472 22.859 40.996 6.610 6.583 1.283 0.447 1.950 – 30.339 10.843 K 0.0066 10.548 3.674 0.753 1.351 0.218 0.217 0.042 0.015 0.064 0.033 – 0.357 L 0.0183 29.516 10.281 2.108 3.781 0.610 0.607 0.118 0.041 0.180 0.092 2.798 – Model Model likelihoods A B C D E F G H I J K L A 0.0006 – 0.348 0.071 0.128 0.021 0.021 0.004 0.001 0.006 0.003 0.095 0.034 B 0.0018 2.871 – 0.205 0.368 0.059 0.059 0.012 0.004 0.017 0.009 0.272 0.097 C 0.0087 14.000 4.876 – 1.793 0.289 0.288 0.056 0.020 0.085 0.044 1.327 0.474 D 0.0049 7.806 2.719 0.558 – 0.161 0.161 0.031 0.011 0.048 0.024 0.740 0.264 E 0.0300 48.419 16.865 3.459 6.202 – 0.996 0.194 0.068 0.295 0.151 4.590 1.640 F 0.0302 48.613 16.933 3.472 6.227 1.004 – 0.195 0.068 0.296 0.152 4.609 1.647 G 0.1547 249.355 86.854 17.811 31.942 5.150 5.129 – 0.348 1.519 0.779 23.639 8.448 H 0.4439 715.903 249.360 51.136 91.707 14.785 14.727 2.871 – 4.362 2.237 67.869 24.255 I 0.1018 164.129 57.169 11.724 21.025 3.390 3.376 0.658 0.229 – 0.513 15.560 5.561 J 0.1985 320.032 111.472 22.859 40.996 6.610 6.583 1.283 0.447 1.950 – 30.339 10.843 K 0.0066 10.548 3.674 0.753 1.351 0.218 0.217 0.042 0.015 0.064 0.033 – 0.357 L 0.0183 29.516 10.281 2.108 3.781 0.610 0.607 0.118 0.041 0.180 0.092 2.798 – The Bayes factors have a row to column interpretation (the value at row i and column j is BFij = MLi/MLj). These values correspond to the analysis with a mutation rate of 5.14% per Myr, after averaging the five smallest acceptance proportions (0.001–0.005%) of the closest simulated data sets to the observed one. The best-fit model is shown in bold. View Large Table 1. Model likelihoods and Bayes factors among the 12 models compared by approximate Bayesian computation Model Model likelihoods A B C D E F G H I J K L A 0.0006 – 0.348 0.071 0.128 0.021 0.021 0.004 0.001 0.006 0.003 0.095 0.034 B 0.0018 2.871 – 0.205 0.368 0.059 0.059 0.012 0.004 0.017 0.009 0.272 0.097 C 0.0087 14.000 4.876 – 1.793 0.289 0.288 0.056 0.020 0.085 0.044 1.327 0.474 D 0.0049 7.806 2.719 0.558 – 0.161 0.161 0.031 0.011 0.048 0.024 0.740 0.264 E 0.0300 48.419 16.865 3.459 6.202 – 0.996 0.194 0.068 0.295 0.151 4.590 1.640 F 0.0302 48.613 16.933 3.472 6.227 1.004 – 0.195 0.068 0.296 0.152 4.609 1.647 G 0.1547 249.355 86.854 17.811 31.942 5.150 5.129 – 0.348 1.519 0.779 23.639 8.448 H 0.4439 715.903 249.360 51.136 91.707 14.785 14.727 2.871 – 4.362 2.237 67.869 24.255 I 0.1018 164.129 57.169 11.724 21.025 3.390 3.376 0.658 0.229 – 0.513 15.560 5.561 J 0.1985 320.032 111.472 22.859 40.996 6.610 6.583 1.283 0.447 1.950 – 30.339 10.843 K 0.0066 10.548 3.674 0.753 1.351 0.218 0.217 0.042 0.015 0.064 0.033 – 0.357 L 0.0183 29.516 10.281 2.108 3.781 0.610 0.607 0.118 0.041 0.180 0.092 2.798 – Model Model likelihoods A B C D E F G H I J K L A 0.0006 – 0.348 0.071 0.128 0.021 0.021 0.004 0.001 0.006 0.003 0.095 0.034 B 0.0018 2.871 – 0.205 0.368 0.059 0.059 0.012 0.004 0.017 0.009 0.272 0.097 C 0.0087 14.000 4.876 – 1.793 0.289 0.288 0.056 0.020 0.085 0.044 1.327 0.474 D 0.0049 7.806 2.719 0.558 – 0.161 0.161 0.031 0.011 0.048 0.024 0.740 0.264 E 0.0300 48.419 16.865 3.459 6.202 – 0.996 0.194 0.068 0.295 0.151 4.590 1.640 F 0.0302 48.613 16.933 3.472 6.227 1.004 – 0.195 0.068 0.296 0.152 4.609 1.647 G 0.1547 249.355 86.854 17.811 31.942 5.150 5.129 – 0.348 1.519 0.779 23.639 8.448 H 0.4439 715.903 249.360 51.136 91.707 14.785 14.727 2.871 – 4.362 2.237 67.869 24.255 I 0.1018 164.129 57.169 11.724 21.025 3.390 3.376 0.658 0.229 – 0.513 15.560 5.561 J 0.1985 320.032 111.472 22.859 40.996 6.610 6.583 1.283 0.447 1.950 – 30.339 10.843 K 0.0066 10.548 3.674 0.753 1.351 0.218 0.217 0.042 0.015 0.064 0.033 – 0.357 L 0.0183 29.516 10.281 2.108 3.781 0.610 0.607 0.118 0.041 0.180 0.092 2.798 – The Bayes factors have a row to column interpretation (the value at row i and column j is BFij = MLi/MLj). These values correspond to the analysis with a mutation rate of 5.14% per Myr, after averaging the five smallest acceptance proportions (0.001–0.005%) of the closest simulated data sets to the observed one. The best-fit model is shown in bold. View Large Four models obtained higher likelihood support than the rest (Figs 2I, 3G–J). In all four models, the most recent populations originated before the LGM in two waves of colonization. These models had BF values in the range of 7.5–360, 4.3–3800, 8.6–1560 and 18.0–105 for mutation rates of 5.14, 10, 20 and 30% per Myr, respectively. These results provide strong support for the presence of at least three waves of sequential colonization (and extinction/replacement) of Crimean red deer by unrelated Western populations, and reject that the present-day population originated from a local pre-LGM population or even a post-LGM population that existed until 2000 yr BP. For mutation rates of 5.14 and 10% per Myr, the BF was > 2.3 for the models with the admixed pre-LGM population (Fig. 3H, J), and it was > 1.5 for the models G and I where the pre-LGM population was panmictic (Fig. 3G, I). However, for mutation rates of 20 and 30% per Myr, the converse occurred (BF = 0.12 and 0.06, respectively). Comparison of models also indicates that the support for a Western origin of post-LGM populations in at least two colonization waves is strong and consistent. Finally, the POD analysis showed low statistical power for several models (Supporting Information, Table S8) ranging from 0.22 to 0.48, which rose to 0.28–0.83 when only BFs equal to or larger than the observed ones were taken into account. The lack of statistical power resulted from ‘competition’ between similar models. An exploratory analysis showed that if similar models (e.g. models of Fig. 3G, H) were mixed into one, the statistical power rose to over 0.9 even for the least supported scenarios. Therefore, the POD analysis showed that the model choice procedure is able to separate correctly between scenarios, although it is less sensitive to distinguishing between very similar models. DISCUSSION Phylogeography of Cervus lineages including Crimean samples We analysed the most sequence-rich set of cytochrome b from different Cervini using the greatest variety of methods so far. The obtained global tree topologies are generally in agreement with results of other authors (Pitra et al., 2004; Hassanin et al., 2012) – Supporting Information, Figure S2. Divergence times for particular lineages of Cervini estimated by us are close to those by Gilbert, Ropiquet & Hassanin (2006), although slightly younger (Supporting Information, Fig. S3). Our phylogenetic study clearly recovered two monophyletic clades of Cervus, corresponding to the eastern and western part of its distribution (Randi et al., 2001; Ludt et al., 2004; Pitra et al., 2004; Skog et al., 2009). The analyses of a large number of deer samples, including fossils from Crimea, also enabled us to infer branching order and migration routes of Cervus in more detail. The significant grouping of the Eastern lineage with Przewalskium from China suggests that this lineage emerged somewhere in Asia. The phylogenies are based on an mtDNA marker, which reflects the history of the maternal lineage. Therefore, nuclear markers should be studied to verify the inferred relationships. About 1 Mya, the Eastern lineage separated into sika (C. nippon) and wapiti (C. canadensis) clades, which further differentiated in Central and Eastern Asia (Fig. 4; Supporting Information, Figs S3, S4). The former colonized Japanese islands (Nagata et al., 1999; Groves, 2005) and the latter reached north-eastern Siberia in the pre-LGM period (Meiri et al., 2014). Wapiti deer carrying mtDNA haplotypes belonging to the E3 clade showed the most widespread distribution because they not only colonized North America but also migrated as far west as Crimea, where they survived to ~35 ka. This implies long-distance migrations of the eastern deer, originally confined into Central and Eastern Asia. The colonization of North America was a single event (Meiri et al., 2014) (Fig. 4) and probably started from north-eastern Siberia (Supporting Information, Fig. S3). Our molecular dating indicates that the migration occurred between 12.4 and 75.9 ka. The deer that crossed the land bridge across the present Bering Strait to North America carried at least two mitochondrial haplotypes, represented by clades E1 and E3. The ancestor of the Western lineage most likely lived in Central Asia, as indicated by the earliest diverged W1 clade (Fig. 5; Supporting Information, Fig. S5). This clade includes Yarkand deer (C. elaphus yarkandensis) and Bactrian deer (C. elaphus bactrianus) from China and Uzbekistan, closely related to C. elaphus hanglu from Indian Kashmir (Lorenzini & Garofalo, 2015). Lorenzini & Garofalo (2015) proposed to raise these three deer to species level as Cervus hanglu. The migration of red deer to Europe happened likely about 1 Mya, as implied by the oldest fossils of C. elaphus (Brugal & Croitor, 2007) – Figure 6A. This is supported by our phylogenies, where an early separated clade with Spanish samples (W2) is basal to the rest of the sequences (Supporting Information, Fig. S5). This is also consistent with our estimated split time between clade W1 and the rest: 1.08 (0.57–1.72) Mya (Supporting Information, Fig. S3). The ancient Spanish samples are dated to post-LGM and do not group with other Spanish sequences, which separated much later (the W9 clade), suggesting that they represent an old lineage that became extinct without descendants. A similar placement of these sequences was obtained in the phylogenetic tree by Meiri et al. (2013). The migration route to Europe probably extended across the Near East and Asia Minor (Fig. 6A) because sequences from Iran, Jordan and Turkey clades are closely related with those from Spain (Supporting Information, Fig. S5). Figure 6. View largeDownload slide Inferred migration routes of red deer in Europe. A, the first settlement of Europe from Central Asia across Near East and Asia Minor as well as migrations to the southern region of Eastern Europe including the Crimean Peninsula. B, further migrations from the latter region to Europe and towards the Ural Mountains. C, post-LGM recolonization of the rest territory of Europe including the Crimean Peninsula from Iberian refugium. The role of the Balkan and Apennine refugia is not shown. Figure 6. View largeDownload slide Inferred migration routes of red deer in Europe. A, the first settlement of Europe from Central Asia across Near East and Asia Minor as well as migrations to the southern region of Eastern Europe including the Crimean Peninsula. B, further migrations from the latter region to Europe and towards the Ural Mountains. C, post-LGM recolonization of the rest territory of Europe including the Crimean Peninsula from Iberian refugium. The role of the Balkan and Apennine refugia is not shown. The next diverged W5 clade includes sequences from the Caucasus, Krasnodarsky Krai, Voronezh and Crimea, which suggests that populations from the Near East also migrated north along the western or/and eastern seashore of the Caucasus and inhabited the southern part of Eastern Europe, including Crimea (Fig. 6A). The first migration may have occurred around 0.6 Mya because the oldest red deer from the Caucasus is dated to 583 ka (Blackwell et al., 2005). This agrees with our calculation of separation time between the lineage including Turkey and Iran sequences and other clades, that is 512 (265–832) ka (Supporting Information, Fig. S3). Samples from Crimea and the Caucasus are separated into different clades with short internal branches, which may indicate a quite diversified ancestral population or several rapid migration events to the north. Some deer from the region migrated to the northeast (Fig. 6A) as indicated by ancient sequences from the Urals, identical to those from the Caucasus (Supporting Information, Fig. S5). At least four westward migrations could have occurred from the southern part of Eastern Europe (Fig. 6B) or different haplotypes were carried by a larger population during a smaller number of migration events. This is suggested by clades W6–8 including pre-LGM sequences from Crimea and the Caucasus and from different regions of Europe (Fig. 5; Supporting Information, Fig. S5). In the W6 and W8 clades, these sequences occupy a basal position to others (Supporting Information, Fig. S5). The separation of these clades appears to have begun roughly 468 (244–758) ka (Supporting Information, Fig. S3). Drastic cooling during the LGM (23–19 ka) caused the distribution of red deer in Europe to be limited to the Balkan, Italian and Iberian Peninsulas, south-western France and the region east of the Carpathians (Banks et al., 2008; Sommer et al., 2008; Skog et al., 2009; Sommer & Zachos, 2009; Meiri et al., 2013). After the LGM with the start of Bølling-Allerød warming (15–13 ka), deer expanded from refugial areas into other regions of Europe. Our results are in agreement with those studies. In clade W9 (corresponding to haplogroup A), the basal position is occupied by Spanish sequences dated from c. 19 ka to the Neolithic (Meiri et al., 2013) – Supporting Information, Figure S5. Other Spanish sequences are younger and are intermingled in the tree with sequences from Western and Central Europe, in line with expectations under a scenario where the red deer that colonized other regions of Europe originated in the Iberian Peninsula (Ludt et al., 2004; Sommer et al., 2008; Skog et al., 2009; Sommer & Zachos, 2009; Niedziałkowska et al., 2011; Meiri et al., 2013) – Figure 6C. Sequences from Crimean deer are separated in the W9 clade, which suggests two independent waves of migration into the Crimean Peninsula or one colonization event by a deer population that carried different haplotypes. A subcluster in clade W6 (corresponding to C haplogroup) with sequences from Southern and Central Europe may correspond to the Balkan refugium. Sequences from Turkey are also present in this clade, which could imply a migration between Europe and Asia Minor. Alternatively, deer from the Balkan refugium could have spread further to the south-east. The W8 clade (corresponding to B haplogroup), besides harbouring sequences from Corsican and Tunisian red deer (C. elaphus corsicanus and C. elaphus barbarus), also includes sequences from Hungary and the Swiss Alps. These most likely represent human translocations from Italy. Detailed population history of Crimean red deer The goal of our study was to investigate the role of the Crimean Peninsula as a possible Eastern European refugium. In agreement with Stankovic et al. (2011), we found that red deer from both the Eastern and Western lineages were present in the Crimean Peninsula before the LGM, even though these populations have otherwise been isolated since their divergence (Kuwayama & Ozawa, 2000; Mahmut et al., 2002). Our ABC analyses, assuming the mutation rates of 5.14 and 10% per Myr, showed that the pre-LGM Crimean population was probably the result of a recent mixture of these lineages. However, for higher mutation rates, the best-supported model (G) assumes that the pre-LGM population was panmictic. This discrepancy results from the model selection procedure, where the POD analysis had difficulties distinguishing between very similar models. Based solely on the ABC results, we cannot therefore be certain which of the scenario regarding the pre-LGM population is more probable. However, the scenario presented in model H is in agreement with other results: the early divergence estimate for the two red deer lineages (c. 1.54 Mya), the distribution of pre-LGM specimens in separate clades, E3 and W6 from the Eastern and Western lineages, and radiocarbon dates and the distribution of red deer samples in cave layers. The specimens from clade W6 were usually found in the lower part of sediments in Emine-Bair-Khosar Cave (Supporting Information, Fig. S1). This layer (I) was dated with pollen to the Eemian stage (MIS5e) when the climate was warmer than now (Avdeienko, 2015) – Supporting Information, Table S2. The oldest deer remains (J46 and J47) with 14C radiocarbon dates close or above the upper limit of this method also belong to clade W6. On the other hand, samples from the eastern population, located in clade E3, were frequent in the upper layer (H), dated to the Middle Pleniglacial (MIS3) by 14C and pollen (Avdeienko, 2015) – Supporting Information, Figure S1 and Table S2. The deer lived in a boreal climate, which was cooler than now, and became extinct before the LGM when the accumulation of allochthonous material in the cave stopped. The youngest pre-LGM sample J28 (28332 cal yr BP) again belongs to clade W6. The data suggest that the deer from clade W6 dominated in Crimea earlier, next the deer from clade E3 appeared and again the deer from clade W6 returned (Fig. 2). However, some radiocarbon dates of samples from these two lineages overlap. Therefore, we cannot exclude that the different deer could interact in some periods (e.g. between 47400 and 43700 cal yr BP), although there were also periods when only one deer clade existed, for example before c. 50000 cal yr BP or between 37300 and 34900 cal yr BP. We did not observe a genetic continuity between pre- and post-LGM samples from Crimea (Fig. 2), which suggests a local extinction hypothesis. This is also supported by the lack of red deer remains from the LGM period in Emine-Bair-Khosar Cave, where such remains are generally abundant (Vremir & Ridush, 2005). However, it cannot be excluded that bone accumulation was impossible in the MIS 2 due to a snow-ice plug that covered the entrance to the cave (Ridush et al., 2013). In this time, the Scandinavian Ice Sheet reached its maximum extent, the climate became colder and drier, the sea level dropped and many temperate species retreated to refugial areas in the southern peninsulas (Hewitt, 2004; Hubberten, 2004; Clark et al., 2009). Environmental conditions in Crimea also deteriorated at the end of MIS 3. Pollen, palaeontological and stable isotope analyses indicated that the climate in the LGM was cold and dry, which reduced broadleaved trees and extended steppes (Gerasimenko, 2007; Markova, 2011; Gąsiorowski et al., 2014). Those unfavourable conditions could have contributed to the extinction of the pre-LGM red deer population. However, the results of Markova (2011) revealed that despite climate deterioration, Crimea was a refugium for small mammals. Therefore, we cannot exclude that climate was the only factor that influenced the red deer extinction. The post-LGM samples are distributed into three clades that belong to the Western lineage and were not descendants of the pre-LGM population, as is indicated by our phylogenetic and ABC analyses. Thus, the pre-LGM population became extinct and was replaced by immigrants. Our data suggest that the population replacement started around 27–14 ka when the pre-LGM deer disappeared from Crimea (Fig. 2). The oldest immigrant population is represented by samples J16 and J27 (c. 13900 cal yr BP), which were probably derived from Russian and Caucasian populations (assigned in Fig. 5 as J16* and in Supporting Information, Fig. S5 as JX966135* in clade W5). They could get into Crimea during the Bølling-Allerød interstadial along the Black Sea seashore from the Transcaucasia refugim, but the lineage vanished in the next Younger Dryas cooling. Younger samples (from Kr to BSN13 dated to c. 12000–2400 cal yr BP, in Fig. 2) belong to clade W7 and are grouped together with sequences from Serbia and Italy and pre- and post-LGM specimens from England. This suggests that Crimea was colonized from the Balkans in the beginning of Holocene. Noticeably, in the late Pleistocene and early Holocene, the Black Sea level was lower than today (Fig. 1). As a consequence, the Balkans and the Caucasus were connected to Crimea by land, which made the migration between these areas possible (Stanko, 2007; Demidenko, 2008). This is also supported by similarities in material culture from the Late Palaeolithic and Mesolithic sites in Bulgaria, Romania and Crimea, suggesting direct contacts among humans from these localities (Stanko, 2007). Despite more stable and favourable climatic conditions since the beginning of the Holocene, the post-LGM red deer population from clade W7 did not survive to modern times in Crimea. This population disappeared c. 2.5 ka and was not an ancestor of contemporary deer. Causes for this Holocene extinction are less clear. Based on soil development and pollen analysis, Cordova & Lehman (2005) reconstructed the environmental changes in Crimea during the Holocene. Their study revealed climate fluctuations and two events of dry and cool phases with a decrease in arboreal vegetation and spreading steppe taxa, dated to 11–7.5 and 4.2–3.5 ka, respectively. This climatic deterioration, however, is unlikely to have caused the extinction of the Holocene population, given that similar conditions are present today in some parts of Crimea, where red deer exist. Cordova & Lehman (2005) also found an increase in phryganoid plants, weeds and other cultivated species around 3.5 ka, suggesting a human impact on the Crimean environment. Between seventh and fifth centuries BC (2.65–2.45 ka), Greeks started to colonize Crimea. One of the earliest settlements on the northern Black Sea coast was established c. 700 yr BC (2.65 ka) on Berezan Island, which was then a part of the mainland (Dolukhanov & Shilik, 2007). The Greek colonization contributed not only to increasing amounts of pastoral plants but also to intensive deforestation, which took place c. 2.5 ka (Cordova & Lehman, 2005; Cordova, 2007). This human activity corresponds well with the observed Holocene red deer extinction, suggesting that intensive landscape management led to red deer habitat fragmentation and possibly to its extinction. Humans have been postulated as a cause of faunal extinction in many other studies (Burney et al., 2004; Burney & Flannery, 2005; Lorenzen et al., 2011; Rawlence et al., 2015). Some of these concentrated on the overkill hypothesis, while others suggested less direct reasons such as habitat alteration, biological invasion or diseases (Burney & Flannery, 2005). The second recolonization of the peninsula appears to have started ~1000 years after the post-LGM extinction of red deer, between the Late Antic and Mediaeval periods, as is supported by two Crimean specimens (SOKIL and PR2) dated to 1.2–0.9 cal ka (eighth to 11th centuries AD) – Figure 2. These, together with present samples, belong to Western haplogroup A (clade W9). This haplogroup originated on the Iberian Peninsula and spread across the whole of Europe (Skog et al., 2009; Niedziałkowska et al., 2011; Meiri et al., 2013). This suggests that Crimea was probably colonized naturally by immigrant red deer from Western Europe. However, translocations by humans cannot be excluded, as it is well known that this happened in Sardinia, Corsica, Crete and Ireland during the Neolithic period (Hajji et al., 2008; Carden et al., 2012; Harris, 2014). The reasons for such management are not well known, but studies by Sykes et al. (2006, 2011) and Davis & MacKinnon (2009) have showed that fallow deer (Dama dama) was introduced to Britain and Portugal by the Romans, probably for keeping them in special parks (vivaria), which enhanced their owner’s status. In the tenth century AD, Crimea was a part of the Eastern Roman Empire. Historical and archaeological studies suggest that deer was a totemic animal for Romans and a symbol of wealth and power. Moreover, hunting for red deer in Roman Britain was an expression of status (Allen et al., 2014). Therefore, a potential extinction of natural red deer population due to Greek activities in Crimea may have induced the need for its reintroduction. More studies on ancient red deer populations, especially from potential migration routes from Iberia to Crimea, are crucial to resolve the question about the origin of the contemporary red deer population in Crimea. Comparison of the 12 models with ABC analyses revealed that after the LGM, red deer likely colonized the Crimean Peninsula in two independent waves: at the end of the Pleistocene and in the Holocene. Those events were related to extinctions and replacements of populations. Similar events have been reported for the collared lemming, where population turnovers were probably caused by climate fluctuations and changes in vegetation (Brace et al., 2012; Palkopoulou et al., 2016). Red deer, however, has a wider environmental tolerance and is generally well adapted to various environments (Straus, 1981; Geist, 1998). Thus, climate as a single factor driving extensive changes in its population may be less likely than for small mammals. Correspondingly, red deer as an important game species has been hunted, translocated and introduced since ancient times (Hartl, Zachos & Nadlinger, 2003; Zachos & Hartl, 2011). We find no support for the hypothesis that Crimea was a Pleistocene refugium for red deer. At least two extinction events in the late Quaternary are in conflict with the classical refugium definition. Considering that our analyses were based on an mtDNA marker and thus reflect only the maternal lineage history, it would be interesting to also analyse nuclear markers to confirm our conclusions on red deer population history in Crimea. Nonetheless, the Crimean samples analysed here provide significant insights into the phylogeography and migration routes of red deer at a global scale. SUPPORTING INFORMATION Additional Supporting Information may be found in the online version of this article at the publisher’s web-site: Table S1. Detailed information about samples used in the study. Radiocarbon dates obtained in previous studies are marked by asterisks. GenBank accession numbers in italics indicate partial cytb sequences published by Stankovic et al. (2011). Table S2. Chronostratigraphy of layers in Emine-Bair-Khosar Cave. Table S3. Nucleotide substitution models for alignment sets used in three programs. In all models, we assumed five discrete categories in the case of gamma distributed rates. Table S4. Summary statistics used in ABC analysis. N – number of individuals; h – number of haplotypes; S – segregating sites; π – nucleotide diversity; PairDiff – pairwise differences between groups. Table S5. Model likelihoods and Bayes factors among the 12 models compared by approximate Bayesian computation. The Bayes factors have a row to column interpretation (the value at row i and column j is BFij = MLi/MLj). These values correspond to the analysis with a mutation rate of 10.0% per million years, after averaging the five smallest acceptance proportions (0.001 to 0.005% of the closest simulated datasets to the observed one). The best-fitted model was bolded. Table S6. Bayes factors among the 12 models compared by approximate Bayesian computation. The Bayes factors have a row to column interpretation (the value at row i and column j is BFij = MLi/MLj). These values correspond to the analysis with a mutation rate of 20.0% per million years, after averaging the ten smallest acceptance proportions (0.001–0.01% of the closest simulated datasets to the observed one). The number of acceptance proportions averaged for this result was higher than in tables with lower mutation rates due to many models presented a likelihood of 0.0 in the first five proportions. The best-fitted model was bolded. Table S7. Bayes factors among the 12 models compared by approximate Bayesian computation. The Bayes factors have a row to column interpretation (the value at row i and column j is BFij = MLi/MLj). These values correspond to the analysis with a mutation rate of 30.0% per million years, after averaging the ten smallest acceptance proportions (0.001 to 0.005% of the closest simulated datasets to the observed one). The number of acceptance proportions averaged for this result was higher than in tables with lower mutation rates due to several models presented a likelihood of 0.0 in the first five proportions. The best-fitted model was bolded. Table S8. Rate of error and statistical power in the model choice analysis by approximate Bayesian computation, as assessed by pseudo observed data sets (PODs). Figure S1. Cross section through layers in Emine-Bair-Khosar Cave. The position of red deer remains was marked. Uncalibrated radiocarbon dates of these samples and membership of DNA from the samples to phylogenetic clades were given. Figure S2. MrBayes tree based on the longer alignment of Cervini members with simplification of Eastern (clade E) and Western (clade W) lineages of red deer. The asterisk indicates that a given sequence is represented by more than one accession number. The shading of circles at nodes corresponds to the number of methods (out of five) which recovered the given node. Numbers at nodes, in the order shown, correspond to: posterior probabilities estimated in MrBayes (MB) and PhyloBayes (PB), bootstrap values obtained in PhyML(PH), TreeFinder (TF) and PAUP by neighbor-joining methods (NJ). Values of the posterior probabilities and bootstrap percentages ≤ 0.50 and 50%, respectively, were omitted or indicated by a dash ‘–’. Figure S3. Simplified maximum clade credibility tree obtained in BEAST for the longer alignment of Cervini members. Median ages expressed in ka BP are shown for all nodes. Sequences obtained from Crimea specimens are bolded. The Eastern (clade E) and Western (clade W) lineages of red deer were indicated. Figure S4. MrBayes tree based on the shorter alignment for the Eastern lineage of red deer. Clades E1-E3 represent C. canadensis (wapiti), whereas E4-E7 C. nippon (sika). Sequences newly obtained in this study from Crimea specimens are bolded. The asterisk indicates that a given sequence is represented by more than one accession number. The shading of circles at nodes corresponds to the number of methods (out of five) which recovered the given node. A probable age or calibrated dates are shown in square brackets. Numbers at nodes, in the order shown, correspond to: posterior probabilities estimated in MrBayes (MB) and PhyloBayes (PB), bootstrap values obtained in PhyML (PH), TreeFinder (TF) and PAUP by neighbor-joining methods (NJ). Values of the posterior probabilities and bootstrap percentages ≤ 0.50 and 50%, respectively, were omitted or indicated by a dash ‘–’. Figure S5. MrBayes tree based on the shorter alignment for the Western lineage of red deer. Sequences newly obtained in this study from Crimea specimens are bolded. The asterisk indicates that a given sequence is represented by more than one accession number. The shading of circles at nodes corresponds to the number of methods (out of five) which recovered the given node. A probable age or calibrated dates are shown in square brackets. Numbers at nodes, in the order shown, correspond to: posterior probabilities estimated in MrBayes (MB) and PhyloBayes (PB), bootstrap values obtained in PhyML (PH), TreeFinder (TF) and PAUP by neighbor joining methods (NJ). Values of the posterior probabilities and bootstrap percentages 0.50 and 50%, respectively, were omitted or indicated by a dash ‘–’. Box S1. Input file for scenario A with mutation rate 5.14% per million years for final simulations. Box S2. Input file for scenario B with mutation rate 5.14% per million years for final simulations. Box S3. Input file for scenario C with mutation rate 5.14% per million years for final simulations. Box S4. Input file for scenario D with mutation rate 5.14% per million years for final simulations. Box S5. Input file for scenario E with mutation rate 5.14% per million years for final simulations. Box S6. Input file for scenario F with mutation rate 5.14% per million years r for final simulations. Box S7. Input file for scenario G with mutation rate 5.14% per million years for final simulations. Box S8. Input file for scenario H with mutation rate 5.14% per million years for final simulations. Box S9. Input file for scenario I with mutation rate 5.14% per million years for final simulations. Box S10. Input file for scenario J with mutation rate 5.14% per million years for final simulations. Box S11. Input file for scenario K with mutation rate 5.14% per million years for final simulations. Box S12. Input file for scenario L with mutation rate 5.14% per million years for final simulations. ACKNOWLEDGEMENTS We are very grateful to three anonymous reviewers for their valuable comments and insightful remarks that have significantly improved the paper. We thank prof. Jan Burdukiewicz who supplied dates for prehistoric periods and archaeological cultures that were ascribed for some ancient samples. We are thankful to Dr R. Karasiewicz-Szypiorski for helpful discussions about Crimean history. The research for this paper was supported by Ministry of Science and Higher Education grant no. N N303 821140 and by the EU through the European Social Fund, contract number UDA-POKL.04.01.01-00-072/09-00. REFERENCES Allen MG , Baker K , Carden R , Madgwick R . 2014 . Chasing Sylvia’s Stag: placing deer in the countryside of Roman Britain . In: Baker K , Carden R , Madgwick R , eds. Deer and people . Oxbow, Oxford, UK: Windgather Press , 174 – 186 . Avdeienko Y . 2015 . The Last Interglacial vegetation and climate in the karstic areas of the Crimea and the Middle Dniester Basin . Physical Geography and Geomorphology (Kyiv) 4 : 102 – 109 . Banks W , Derrico F , Peterson A , Kageyama M , Colombeau G . 2008 . Reconstructing ecological niches and geographic distributions of caribou (Rangifer tarandus) and red deer (Cervus elaphus) during the Last Glacial Maximum . Quaternary Science Reviews 27 : 2568 – 2575 . Google Scholar CrossRef Search ADS Beaumont MA , Zhang W , Balding DJ . 2002 . Approximate Bayesian computation in population genetics . Genetics 162 : 2025 – 2035 . Bennett K , Provan J . 2008 . What do we mean by ‘refugia’ ? Quaternary Science Reviews 27 : 2449 – 2455 . Google Scholar CrossRef Search ADS Blackwell BA , Liang S , Golovanova LV , Doronichev VB , Skinner AR , Blickstein JI . 2005 . ESR at Treugol’naya Cave, Northern Caucasus Mt., Russia: dating Russia’s oldest archaeological site and paleoclimatic change in Oxygen Isotope Stage 11 . Applied Radiation and Isotopes 62 : 237 – 245 . Google Scholar CrossRef Search ADS Brace S , Palkopoulou E , Dalén L , Lister AM , Miller R , Otte M , Germonpré M , Blockley SPE , Stewart JR , Barnes I . 2012 . Serial population extinctions in a small mammal indicate Late Pleistocene ecosystem instability . Proceedings of the National Academy of Sciences of the United States of America 109 : 20532 – 20536 . Google Scholar CrossRef Search ADS Brugal JP , Croitor R . 2007 . Evolution, ecology and biochronology of herbivore associations in Europe during the last 3 million years . Quaternaire 18: 129 – 152 . Burney DA , Burney LP , Godfrey LR , Jungers WL , Goodman SM , Wright HT , Jull AJ . 2004 . A chronology for late prehistoric Madagascar . Journal of Human Evolution 47 : 25 – 63 . Google Scholar CrossRef Search ADS Burney DA , Flannery TF . 2005 . Fifty millennia of catastrophic extinctions after human contact . Trends in Ecology & Evolution 20 : 395 – 401 . Google Scholar CrossRef Search ADS Cameron RAD , Pokryszko BM , Horsák M . 2013 . Forest snail faunas from Crimea (Ukraine), an isolated and incomplete Pleistocene refugium . Biological Journal of the Linnean Society 109 : 424 – 433 . Google Scholar CrossRef Search ADS Carden RF , McDevitt AD , Zachos FE , Woodman PC , O’Toole P , Rose H , Monaghan NT , Campana MG , Bradley DG , Edwards CJ . 2012 . Phylogeographic, ancient DNA, fossil and morphometric analyses reveal ancient and modern introductions of a large mammal: the complex case of red deer (Cervus elaphus) in Ireland . Quaternary Science Reviews 42 : 74 – 84 . Google Scholar CrossRef Search ADS Clark PU , Dyke AS , Shakun JD , Carlson AE , Clark J , Wohlfarth B , Mitrovica JX , Hostetler SW , McCabe AM . 2009 . The Last Glacial Maximum . Science 325 : 710 – 714 . Google Scholar CrossRef Search ADS Cooper A , Turney C , Hughen KA , Barry W , Mcdonald HG , Bradshaw CJA . 2015 . Abrupt warming events drove Late Pleistocene Holarctic megafaunal turnover . Science Express 349 : 1 – 8 . Cordova CE . 2007 . Holocene Mediterranization of the southern Crimean vegetation: paleoecological records, regional climate change, and possible non-climatic influences . In: Yanko-Hombach V , Gilbert AS , Panin N , Dolukhanov PM , eds. The Black Sea flood question: changes in coastline, climate, and human settlement . Netherlands: Springer , 319 – 344 . Google Scholar CrossRef Search ADS Cordova CE , Lehman PH . 2005 . Holocene environmental change in southwestern Crimea (Ukraine) in pollen and soil records . The Holocene 15 : 263 – 277 . Google Scholar CrossRef Search ADS Criscuolo A . 2011 . morePhyML: improving the phylogenetic tree space exploration with PhyML 3 . Molecular Phylogenetics and Evolution 61 : 944 – 948 . Google Scholar CrossRef Search ADS Darriba D , Taboada GL , Doallo R , Posada D . 2012 . jModelTest 2: more models, new heuristics and parallel computing . Nature Methods 9 : 772 . Google Scholar CrossRef Search ADS Davis S , MacKinnon M . 2009 . Did the Romans bring fallow deer to Portugal ? Environmental Archaeology 14 : 15 – 26 . Google Scholar CrossRef Search ADS Demidenko YE . 2008 . The Early and Mid-Upper Palaeolithic of the North Black Sea region: an overview . Quartär 55 : 99 – 114 . Dolukhanov PM , Shilik KK . 2007 . Environment, sea-level changes, and human migrations in the northern Pontic area during late Pleistocene and Holocene times . In: Yanko-Hombach V , Gilbert AS , Panin N , Dolukhanov PM , eds. The Black Sea flood question: changes in coastline, climate, and human settlement . Netherlands: Springer , 297 – 318 . Google Scholar CrossRef Search ADS Drummond AJ , Suchard MA , Xie D , Rambaut A . 2012 . Bayesian phylogenetics with BEAUti and the BEAST 1.7 . Molecular Biology and Evolution 29 : 1969 – 1973 . Google Scholar CrossRef Search ADS Fagundes NJR , Ray N , Beaumont M , Neuenschwander S , Salzano FM , Bonatto SL , Excoffier L . 2007 . Statistical evaluation of alternative models of human evolution . Proceedings of the National Academy of Sciences of the United States of America 104 : 17614 – 17619 . Google Scholar CrossRef Search ADS Fernández-García JL , Carranza J , Martínez JG , Randi E . 2013 . Mitochondrial D-loop phylogeny signals two native Iberian red deer (Cervus elaphus) lineages genetically different to Western and Eastern European red deer and infers human-mediated translocations . Biodiversity and Conservation 23 : 537 – 554 . Google Scholar CrossRef Search ADS Gąsiorowski M , Hercman H , Ridush B , Stefaniak K . 2014 . Environment and climate of the Crimean Mountains during the Late Pleistocene inferred from stable isotope analysis of red deer (Cervus elaphus) bones from the Emine-Bair-Khosar Cave . Quaternary International 326–327: 243 – 249 . Geist V . 1998 . Deer of the world: their evolution, behaviour, and ecology . Mechanicsburg, PA: Stackpole Books . Gerasimenko N . 2007 . Environmental changes in the Crimean Mountains during the Last Interglacial–middle pleniglacial as recorded by pollen and lithopedology . Quaternary International 164–165 : 207 – 220 . Google Scholar CrossRef Search ADS Gilbert C , Ropiquet A , Hassanin A . 2006 . Mitochondrial and nuclear phylogenies of Cervidae (Mammalia, Ruminantia): systematics, morphology, and biogeography . Molecular Phylogenetics and Evolution 40 : 101 – 117 . Google Scholar CrossRef Search ADS Groves C . 2005 . The genus Cervus in eastern Eurasia . European Journal of Wildlife Research 52 : 14 – 22 . Google Scholar CrossRef Search ADS Guindon S , Dufayard JF , Lefort V , Anisimova M , Hordijk W , Gascuel O . 2010 . New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0 . Systematic Biology 59 : 307 – 321 . Google Scholar CrossRef Search ADS Hajji GM , Charfi-Cheikrouha F , Lorenzini R , Vigne JD , Hartl GB , Zachos FE . 2008 . Phylogeography and founder effect of the endangered Corsican red deer (Cervus elaphus corsicanus) . Biodiversity and Conservation 17 : 659 – 673 . Google Scholar CrossRef Search ADS Hall TA . 1999 . BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT . Nucleic Acids Symposium Series 41: 95 – 98 . Harris K . 2014 . Hunting, performance and incorporation: human–deer encounter in Late Bronze Age Crete . In: Baker K , Carden R , Madgwick R , eds. Deer and people . Oxbow, Oxford, UK: Windgather Press , 48 – 58 . Hartl GB , Zachos F , Nadlinger K . 2003 . Genetic diversity in European red deer (Cervus elaphus L.): anthropogenic influences on natural populations . Comptes Rendus Biologies 326 ( Suppl 1 ): S37 – S42 . Google Scholar CrossRef Search ADS Hassanin A , Delsuc F , Ropiquet A , Hammer C , Jansen van Vuuren B , Matthee C , Ruiz-Garcia M , Catzeflis F , Areskoug V , Nguyen TT , Couloux A . 2012 . Pattern and timing of diversification of Cetartiodactyla (Mammalia, Laurasiatheria), as revealed by a comprehensive analysis of mitochondrial genomes . Comptes Rendus Biologies 335 : 32 – 50 . Google Scholar CrossRef Search ADS Hewitt GM . 1996 . Some genetic consequences of ice ages, and their role, in divergence and speciation . Biological Journal of the Linnean Society 58 : 247 – 276 . Google Scholar CrossRef Search ADS Hewitt GM . 2004 . Genetic consequences of climatic oscillations in the Quaternary . Philosophical Transactions of the Royal Society of London B: Biological Sciences 359 : 183 – 195 ; discussion 195. Google Scholar CrossRef Search ADS Ho SY . 2014 . The changing face of the molecular evolutionary clock . Trends in Ecology & Evolution 29 : 496 – 503 . Google Scholar CrossRef Search ADS Hubberten H . 2004 . The periglacial climate and environment in northern Eurasia during the Last Glaciation . Quaternary Science Reviews 23 : 1333 – 1357 . Google Scholar CrossRef Search ADS Jobb G , von Haeseler A , Strimmer K . 2004 . TREEFINDER: a powerful graphical analysis environment for molecular phylogenetics . BMC Evolutionary Biology 4 : 18 . Google Scholar CrossRef Search ADS Krojerová-Prokešová J , Barančeková M , Koubek P . 2015 . Admixture of eastern and western European Red deer lineages as a result of postglacial recolonization of the Czech Republic (Central Europe) . The Journal of Heredity 106 : 375 – 385 . Google Scholar CrossRef Search ADS Kuwayama R , Ozawa T . 2000 . Phylogenetic relationships among European red deer, wapiti, and sika deer inferred from mitochondrial DNA sequences . Molecular Phylogenetics and Evolution 15 : 115 – 123 . Google Scholar CrossRef Search ADS Kuznetsova MV , Danilkin AA , Kholodova MV . 2012 . Phylogeography of red deer (Cervus elaphus): analysis of mtDNA cytochrome b polymorphism . Biology Bulletin 39 : 323 – 330 . Google Scholar CrossRef Search ADS Kuznetsova M , Volokh A , Domnich V , Tyshkevitch V , Danilkin A . 2007 . Molecular diversity of the red deer, Cervus elaphus (Cervidae), inhabiting East Europe . Vestnik Zoologii 41 : 505 – 509 . [in Russian] Lartillot N , Lepage T , Blanquart S . 2009 . PhyloBayes 3: a Bayesian software package for phylogenetic reconstruction and molecular dating . Bioinformatics 25 : 2286 – 2288 . Google Scholar CrossRef Search ADS Lorenzen ED , Nogués-Bravo D , Orlando L , Weinstock J , Binladen J , Marske KA , Ugan A , Borregaard MK , Gilbert MTP , Nielsen R , Ho SYW , Goebel T , Graf KE , Byers D , Stenderup JT , Rasmussen M , Campos PF , Leonard JA , Koepfli KP , Froese D , Zazula G , Stafford TW , Aaris-Sørensen K , Batra P , Haywood AM , Singarayer JS , Valdes PJ , Boeskorov G , Burns JA , Davydov SP , Haile J , Jenkins DL , Kosintsev P , Kuznetsova T , Lai X , Martin LD , McDonald HG , Mol D , Meldgaard M , Munch K , Stephan E , Sablin M , Sommer RS , Sipko T , Scott E , Suchard MA , Tikhonov A , Willerslev R , Wayne RK , Cooper A , Hofreiter M , Sher A , Shapiro B , Rahbek C , Willerslev E . 2011 . Species-specific responses of Late Quaternary megafauna to climate and humans . Nature 479 : 359 – 364 . Google Scholar CrossRef Search ADS Lorenzini R , Garofalo L . 2015 . Insights into the evolutionary history of Cervus (Cervidae, tribe Cervini) based on Bayesian analysis of mitochondrial marker sequences, with first indications for a new species . Journal of Zoological Systematics and Evolutionary Research 53 : 340 – 349 . Google Scholar CrossRef Search ADS Ludt CJ , Schroeder W , Rottmann O , Kuehn R . 2004 . Mitochondrial DNA phylogeography of red deer (Cervus elaphus) . Molecular Phylogenetics and Evolution 31 : 1064 – 1083 . Google Scholar CrossRef Search ADS Mahmut H , Masuda R , Onuma M , Takahashi M , Nagata J , Suzuki M , Ohtaishi N . 2002 . Molecular phylogeography of the red deer (Cervus elaphus) populations in Xinjiang of China: comparison with other Asian, European, and North American populations . Zoological Science 19 : 485 – 495 . Google Scholar CrossRef Search ADS Markov GG , Kuznetsova MV , Danilkin AA , Kholodova MV . 2012 . Analysis of genetic diversity of red deer (Cervus elaphus L.) in Bulgaria: implications for population conservation and sustainable management . Acta Zoologica Bulgarica 64 : 389 – 396 . Markov GG , Kuznetsova MV , Danilkin AA , Kholodova MV . 2015 . Genetic diversity of the red deer (Cervus elphus L.) in Hungary revealed by cytochrome b gene . Acta Zoologica Bulgarica 67 : 11 – 17 . Markova AK . 2011 . Small mammals from Palaeolithic sites of the Crimea . Quaternary International 231 : 22 – 27 . Google Scholar CrossRef Search ADS Markova AK , Simakova AN , Puzachenko AY . 2009 . Ecosystems of Eastern Europe at the time of maximum cooling of the Valdai glaciation (24–18 kyr BP) inferred from data on plant communities and mammal assemblages . Quaternary International 201 : 53 – 59 . Google Scholar CrossRef Search ADS Meiri M , Lister AM , Collins MJ , Tuross N , Goebel T , Blockley S , Zazula GD , van Doorn N , Dale Guthrie R , Boeskorov GG , Baryshnikov GF , Sher A , Barnes I . 2014 . Faunal record identifies Bering isthmus conditions as constraint to end-Pleistocene migration to the New World . Proceedings of the Royal Society of London B: Biological Sciences 281 : 20132167 . Google Scholar CrossRef Search ADS Meiri M , Lister AM , Higham TF , Stewart JR , Straus LG , Obermaier H , González Morales MR , Marín-Arroyo AB , Barnes I . 2013 . Late-glacial recolonization and phylogeography of European red deer (Cervus elaphus L.) . Molecular Ecology 22 : 4711 – 4722 . Google Scholar CrossRef Search ADS Nagata J , Masuda R , Tamate HB , Hamasaki Si , Ochiai K , Asada M , Tatsuzawa S , Suda K , Tado H , Yoshida MC . 1999 . Two genetically distinct lineages of the sika deer, Cervus nippon, in Japanese islands: comparison of mitochondrial D-loop region sequences . Molecular Phylogenetics and Evolution 13 : 511 – 519 . Google Scholar CrossRef Search ADS Niedziałkowska M , Jędrzejewska B , Honnen AC , Otto T , Sidorovich VE , Perzanowski K , Skog A , Hartl GB , Borowik T , Bunevich AN , Lang J , Zachos FE . 2011 . Molecular biogeography of red deer Cervus elaphus from eastern Europe: insights from mitochondrial DNA sequences . Acta Theriologica 56 : 1 – 12 . Google Scholar CrossRef Search ADS Palkopoulou E , Baca M , Abramson NI , Sablin M , Socha P , Nadachowski A , Prost S , Germonpré M , Kosintsev P , Smirnov NG , Vartanyan S , Ponomarev D , Nyström J , Nikolskiy P , Jass CN , Litvinov YN , Kalthoff DC , Grigoriev S , Fadeeva T , Douka A , Higham TF , Ersmark E , Pitulko V , Pavlova E , Stewart JR , Węgleński P , Stankovic A , Dalén L . 2016 . Synchronous genetic turnovers across Western Eurasia in Late Pleistocene collared lemmings . Global Change Biology 22 : 1710 – 1721 . doi:10.1111/gcb.13214 . Google Scholar CrossRef Search ADS Pérez-Espona S , Pérez-Barbería FJ , Goodall-Copestake WP , Jiggins CD , Gordon IJ , Pemberton JM . 2009 . Genetic diversity and population structure of Scottish Highland red deer (Cervus elaphus) populations: a mitochondrial survey . Heredity 102 : 199 – 210 . Google Scholar CrossRef Search ADS Pitra C , Fickel J , Meijaard E , Groves PC . 2004 . Evolution and phylogeny of old world deer . Molecular Phylogenetics and Evolution 33 : 880 – 895 . Google Scholar CrossRef Search ADS Polziehn RO , Strobeck C . 2002 . A phylogenetic comparison of red deer and wapiti using mitochondrial DNA . Molecular Phylogenetics and Evolution 22 : 342 – 356 . Google Scholar CrossRef Search ADS R Core Team . 2014 . R: a language and environment for statistical computing . Vienna : R Foundation for Statistical Computing . Rambaut A , Suchard M , Xie D , Drummond A . 2014 . Tracer v1. 6 . Ramsey CB , Scott EM , van der Plicht J . 2013 . Calibration for archaeological and environmental terrestrial samples in the time range 26–50 ka cal bp . Radiocarbon 55 : 2021 – 2027 . Google Scholar CrossRef Search ADS Randi E , Mucci N , Claro-Hergueta F , Bonnet A , Douzery EJP . 2001 . A mitochondrial DNA control region phylogeny of the Cervinae: speciation in Cervus and implications for conservation . Animal Conservation 4 : 1 – 11 . Google Scholar CrossRef Search ADS Rawlence NJ , Perry GLW , Smith IWG , Scofield RP , Tennyson AJD , Matisoo-Smith EA , Boessenkool S , Austin JJ , Waters JM . 2015 . Radiocarbon-dating and ancient DNA reveal rapid replacement of extinct prehistoric penguins . Quaternary Science Reviews 112 : 59 – 65 . Google Scholar CrossRef Search ADS Reimer PJ , Bard E , Bayliss A , Beck JW , Blackwell PG , Bronk Ramsey C , Buck CE , Cheng H , Edwards RL , Friedrich M . 2013 . IntCal13 and Marine13 radiocarbon age calibration curves 0–50,000 years cal BP . Radiocarbon 55 : 1869 – 1887 . Google Scholar CrossRef Search ADS Richter J . 2005 . Hasty foragers: the Crimea Island and Europe during the last interglacial . In: Chabai V , Richter J , Uthmeier T , eds. Kabazi II: Last Interglacial occupation, environment and subsistence, Palaeolithic sites of Crimea . Shlyakh: Simferopol-Cologne , 275 – 286 . Ridush B , Stefaniak K , Socha P , Proskurnyak Y , Marciszak A , Vremir M , Nadachowski A . 2013 . Emine-Bair-Khosar Cave in the Crimea, a huge bone accumulation of Late Pleistocene fauna . Quaternary International 284 : 151 – 160 . Google Scholar CrossRef Search ADS Ronquist F , Teslenko M , van der Mark P , Ayres DL , Darling A , Höhna S , Larget B , Liu L , Suchard MA , Huelsenbeck JP . 2012 . MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space . Systematic Biology 61 : 539 – 542 . Google Scholar CrossRef Search ADS Rosvold J , Røed KH , Hufthammer AK , Andersen R , Stenøien HK . 2012 . Reconstructing the history of a fragmented and heavily exploited red deer population using ancient and contemporary DNA . BMC Evolutionary Biology 12 : 191 . Google Scholar CrossRef Search ADS Sanderson MJ . 2002 . Estimating absolute rates of molecular evolution and divergence times: a penalized likelihood approach . Molecular Biology and Evolution 19 : 101 – 109 . Google Scholar CrossRef Search ADS Sandom C , Faurby S , Sandel B , Svenning JC . 2014 . Global late Quaternary megafauna extinctions linked to humans, not climate change . Proceedings of the Royal Society of London B: Biological Sciences 281 : 20133254 . Google Scholar CrossRef Search ADS Sandoval-Castellanos E , Palkopoulou E , Dalén L . 2014 . Back to BaySICS: a user-friendly program for Bayesian Statistical Inference from Coalescent Simulations . PLoS ONE 9 : e98011 . Google Scholar CrossRef Search ADS Skog A , Zachos FE , Rueness EK , Feulner PGD , Mysterud A , Langvatn R , Lorenzini R , Hmwe SS , Lehoczky I , Hartl GB , Stenseth NC , Jakobsen KS . 2009 . Phylogeography of red deer (Cervus elaphus) in Europe . Journal of Biogeography 36 : 66 – 77 . Google Scholar CrossRef Search ADS Sommer RS , Zachos FE . 2009 . Fossil evidence and phylogeography of temperate species: ‘glacial refugia’ and post-glacial recolonization . Journal of Biogeography 36 : 2013 – 2020 . Google Scholar CrossRef Search ADS Sommer RS , Zachos FE , Street M , Jöris O , Skog A , Benecke N . 2008 . Late Quaternary distribution dynamics and phylogeography of the red deer (Cervus elaphus) in Europe . Quaternary Science Reviews 27 : 714 – 733 . Google Scholar CrossRef Search ADS Stanko VN . 2007 . Fluctuations in the level of the Black Sea and Mesolithic settlement of the northern Pontic area . In: Yanko-Hombach V , Gilbert AS , Panin N , Dolukhanov PM , eds. The Black Sea flood question: changes in coastline, climate, and human settlement . Netherlands: Springer , 371 – 385 . Google Scholar CrossRef Search ADS Stankovic A , Doan K , Mackiewicz P , Ridush B , Baca M , Gromadka R , Socha P , Weglenski P , Nadachowski A , Stefaniak K . 2011 . First ancient DNA sequences of the Late Pleistocene red deer (Cervus elaphus) from the Crimea, Ukraine . Quaternary International 245 : 262 – 267 . Google Scholar CrossRef Search ADS Stewart JR , Lister AM , Barnes I , Dalén L . 2010 . Refugia revisited: individualistic responses of species in space and time . Proceedings of the Royal Society of London B: Biological Sciences 277 : 661 – 71 . Google Scholar CrossRef Search ADS Stiller M , Knapp M , Stenzel U , Hofreiter M , Meyer M . 2009 . Direct multiplex sequencing (DMPS) – a novel method for targeted high-throughput sequencing of ancient and highly degraded DNA . Genome Research 19 : 1843 – 1848 . Google Scholar CrossRef Search ADS Straus LG . 1981 . On the habitat and diet of Cervus elaphus . Munibe 33 : 175 – 182 . Swofford DL . 1998 . Phylogenetic analysis using parsimony (* and other methods), Version 4 . Sykes NJ , Baker KH , Carden RF , Higham TFG , Hoelzel AR , Stevens RE . 2011 . New evidence for the establishment and management of the European fallow deer (Dama dama dama) in Roman Britain . Journal of Archaeological Science 38 : 156 – 165 . Google Scholar CrossRef Search ADS Sykes N , White J , Hayes T , Palmer M . 2006 . Tracking animals using strontium isotopes in teeth: the role of fallow deer (Dama dama) in Roman Britain . 80 : 948 – 960 . Taberlet P , Fumagalli L , Wust-Saucy AG , Cosson JF . 1998 . Comparative phylogeography and postglacial colonization routes in Europe . Molecular Ecology 7 : 453 – 464 . Google Scholar CrossRef Search ADS Tamura K , Stecher G , Peterson D , Filipski A , Kumar S . 2013 . MEGA6: Molecular Evolutionary Genetics Analysis version 6.0 . Molecular Biology and Evolution 30 : 2725 – 2729 . Google Scholar CrossRef Search ADS Vremir MM , Ridush B . 2005 . The Emine-Bair-Khosar ‘Mega Trap’ (Ukraine) . Mitteilungen der Kommission für Quartärforschung Österreichischen Akademie der Wissenschaften 14: 235 – 239 . Zachos FE , Hartl GB . 2011 . Phylogeography, population genetics and conservation of the European red deer Cervus elaphus . Mammal Review 41 : 138 – 150 . Google Scholar CrossRef Search ADS © 2017 The Linnean Society of London, Zoological Journal of the Linnean Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

Journal

Zoological Journal of the Linnean SocietyOxford University Press

Published: Oct 19, 2017

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off