Background

Characterized by considerable climatic instability, the Pleistocene was a time of changing ecological conditions that affected both plant and animal life. Glacial-interglacial cycles were irregular and created contrasting landscapes and environments in many regions over time (Hofreiter and Stewart, 2009). These regional and global changes had a profound impact on animal populations, resulting in lineage diversification, range shifts, demographic bottlenecks, and extinctions. Until recently such processes have mainly been traced through the fossil record at a morphological level. However, in the last couple of decades, the analysis of ancient DNA (aDNA) has been adding valuable information, deciphering ancient population dynamics even in the distant past (Lorenzen et al., 2011).

So far, a number of extant as well as extinct animal species have been investigated using the combined tools of morphology and aDNA analysis. The results have in many instances caused revisions of long-standing theories, giving a more complex picture of how animal populations have reacted to climate change (de Bruyn et al., 2011; Hofreiter and Stewart, 2009). Generally, the model of contraction and expansion of population ranges following the glacial cycles holds true for many species (Bennett and Provan, 2008). However, contrary to the previous belief that animals tracked their habitats during range contractions, it has been shown that some populations died out completely as their environment changed (Campos et al., 2010a; Campos et al., 2010b; Dalen et al., 2007). Furthermore, species often seem to have responded in an individualistic manner when the distribution of habitats varied over time, depending on their ecological requirements and relationships with other species (Lorenzen et al., 2011; Shapiro et al., 2004; Stewart, 2008; Stewart et al., 2010).

The unstable conditions of the Pleistocene seem to have proved fatal for many species. Even though the direct causes are still not fully understood, the Late Quaternary witnessed a major extinction of a great number of large terrestrial mammals (Martin, 1999). The extinctions were not synchronized globally, but by the middle of the Holocene a majority of the so-called megafauna (weighing over 44 kg) outside Africa and Southern Asia had disappeared entirely (Stuart, 1999). Whether these extinctions were caused by climate change, human hunting, or a combination of the two, the result was a dramatic reduction in diversity, in both populations and entire species (Barnosky et al., 2004; Hofreiter, 2007). Some animals, like the woolly mammoth (Mammuthus primigenius) and the giant deer (Megaloceros giganteus), did survive in remote regions for longer than previously supposed (Guthrie, 2004; Stuart et al., 2004; Vartanyan et al., 1993), but they too eventually faced extinction. Another member of this lost fauna, often depicted in Palaeolithic art, was the cave lion (Guthrie, 2005).

Both morphological and mitochondrial DNA (mtDNA) data suggest that the cave lion was a sister species to the African lion (Burger et al., 2004; Mazak, 2010; Sotnikova and Nikolskiy, 2006). However, its relationship to living lions is still not fully resolved and the designation of sub-species status is common (Burger et al., 2004). We have chosen to treat the cave lion as a separate species (Panthera spelaea) following the recent study by Stuart & Lister (Stuart and Lister, 2011).

Compared to extant African lions, cave lions were on average larger (Turner, 1997), with notable differences in cranial and dental morphology (Mazak, 2010; Sotnikova and Nikolskiy, 2006). In other respects their bodily form seems to have been similar to the modern species (Christiansen, 2008). The cave lion’s main prey species probably included horse (Equus ferus), steppe bison (Bison priscus), reindeer (Rangifer tarandus), giant deer (Megaloceros giganteus), red deer (Cervus elaphus), and musk ox (Ovibos moschatus) (Stuart and Lister, 2011). Fossil evidence from Alaska has confirmed bison among its prey (Guthrie, 1990). Stable isotope analyses of cave lion remains from Europe indicate a diverse diet, but also specialization in some individuals on reindeer or juvenile cave bears (Bocherens et al., 2011).

Despite its vast geographical distribution, the cave lion finally went extinct around 14,000 years before present (cal a BP) (Stuart and Lister, 2011). The youngest find in Europe comes from southern Germany and is dated to 14,378 ± 756 cal a BP (Stuart and Lister, 2011). In Siberia, a similar age has been found for a specimen from the Lena delta 14,640 ± 714 cal a BP (Stuart and Lister, 2011), and the latest available date on cave lion, 13,290 ± 549 cal a BP, comes from Fairbanks Creek, Alaska (Stuart and Lister, 2011).

The first aDNA analysis of the cave lion showed it to be closely related to modern African and Asiatic lions, yet genetically distinct (Burger et al., 2004). This was further confirmed by Barnett et al. (2009) with a more extensive dataset from across the species’ range. The latter study also described the American lion (Panthera atrox) as a separate clade and revealed an apparent reduction in genetic diversity across most of the range sometime after 52,000 cal a BP. Nine haplotypes were observed in twelve specimens before that time, in contrast to four among 18 younger specimens. Since the latter four haplotypes were all closely related, with one (haplotype B) at very high frequency, a widespread demographic bottleneck was hypothesized. This scenario is in agreement with a recent large-scale study on radiocarbon-dated cave lion specimens where several temporal gaps in the fossil record were observed, which could indicate local extinctions or declines in population size (Stuart and Lister, 2011).

The cave lion is uniquely well recorded in terms of dated remains among the top predators of the Late Pleistocene (Stuart and Lister, 2011). This consequently offers an insight into the ecology of the Late Pleistocene fauna, unparalleled by data from any other large carnivore. To further investigate the past population dynamics in the cave lion, we analysed mitochondrial DNA variation from a total of 48 specimens collected throughout the species’ former distribution. By including more locations and dates and combining population genetics and phylogenetic analysis with powerful computational analyses, including approximate Bayesian computation (ABC) from serial coalescent simulations and Markov Chain Monte Carlo (MCMC), this study aims to obtain a better inference of the nature of demographic change during the Late Pleistocene both in terms of timing and extent.

Methods

Sampling and DNA extraction

Samples of cave lion (n = 29) from a range of locations across Eurasia and North America were kindly provided by colleagues at various institutions (figure 1, table S2 and Acknowledgements). The samples consisted of bone and tooth fragments, ranging from solid pieces to coarse powder, and varying in quality and amount. All samples had already been radiocarbon dated at the Oxford Radiocarbon Accelerator Unit (ORAU), with ages spanning between >62,000 and 14,000 cal a BP (Stuart and Lister, 2011).

Figure 1 

Geographic distribution of the Pleistocene cave lion. Origins of successfully amplified samples from this study (black circles), and from the study by Barnett et al. (2009) (white circles). The three geographic regions are highlighted; Europe (blue), Siberia (red) and Alaska plus Yukon (yellow).

DNA extractions were performed in a designated ancient DNA laboratory at the Swedish Museum of Natural History in Stockholm. Using a drill (Dremel) or mortar and pestle, about 50 mg of bone from each sample was ground to a fine powder. DNA was extracted using a silica-based method. Roughly 15–50 mg bone powder from each sample was incubated overnight under motion at 55°C in 715 µl extraction buffer (0.45M EDTA, 0.1M UREA, 150 µg proteinase K). Following this, the samples were centrifuged at 2,300 rpm for 5 minutes and the supernatants were collected and concentrated using 30K MWCO Vivaspin filters (Sartorius). The concentrated supernatants were subsequently purified and eluted following the procedure in Yang et al. (Yang et al., 1998).

Amplification and Sequencing

Following Barnett et al. (2009), two regions of the mitochondrial genome were amplified: a 143 base pairs (bp) long fragment of ATP8 and a ~215 bp long fragment of the control region (CR). The latter consisted of two combined and partially overlapping amplicons (without primers); ~101 and ~126 bp long, referred to as CRA and CRB respectively. Primer sequences are given in table S3.

Polymerase chain reactions (PCRs) were prepared in 25 µl volumes, using 1mM MgCl2 (Qiagen), 0.2 mM dNTPs, 0.2 µM of each primer, 1X PCR-buffer (Qiagen), 0.1mg/ml BSA, 2 units of Hotstar Taq (Qiagen) and 2 µl of DNA extract. Amplifications started with a 10 min denaturation step at 95°C, followed by 55 cycles of denaturation at 94°C for 30 s, annealing for 30 s at 56°C for ATP8 and at 46°C for CRA and CRB, followed by extension at 72°C for 30 s. A final extension step at 72°C for 7 min was included at the end of the procedure. Gel electrophoresis was used to detect and confirm amplification success, by applying 5 µl PCR product on a 1.5% agarose gel stained with ethidium bromide. The gel was run in 0.5X TBE buffer and visualized under UV-light. The PCR-products were then treated with EXO-SAP enzymes (1:4) to remove excess dNTPs and unincorporated primers.

Sequencing reactions on both the heavy and light strands were performed using the BigDye Terminator kit ver.1.1 (Applied Biosystems Inc.), followed by purification with the DyeEx 96 Kit (Qiagen). Analysis was then performed on an ABI 3130xl Genetic analyzer (Applied Biosystems Inc.). The resulting sequences were aligned using the Staden Package (Staden, 1999), and consensus sequences verified by eye.

Amplification products that were difficult to sequence, for example those that appeared as double or triple bands after electrophoresis, were cloned with a TOPO TA Cloning Kit (Invitrogen). Cloning was also made on products where direct sequencing resulted in sequences of nuclear mitochondrial inserts (numts), widely occurring among felids (Kim et al., 2006) and with known sequence variants previously reported for cave lions (Barnett et al., 2009). Colonies were grown on agar plates with ampicilin overnight at 37°C, and transferred to distilled water. PCR and sequencing of eight clones per original PCR product were subsequently performed as described above.

Precautions and Authenticity

Special precautions were taken in all steps during the extraction to avoid contamination. The pre-PCR preparation and extraction of the bone samples was performed in a separate facility for aDNA work, with high standards of sterility. All tools and surfaces were carefully washed with 0.5% sodium hypochlorite as well as UV-irradiated, and the equipment was treated in this way between handling of the different samples. Bone surfaces were removed before drilling and contact with the samples was minimised. Extraction series were kept small in order to avoid cross-contamination, and for every fourth sample one extraction blank was included. Furthermore, four samples were extracted twice and one PCR blank was included for every sixth sample to detect possible contamination. At least two PCR products were sequenced in order to assess sequence authenticity and filter for nucleotide misincorporations resulting from post-mortem DNA damage (Hofreiter et al., 2001).

Data analysis

The successful consensus sequences were combined with the dataset from Barnett et al. (2009) to obtain a total of 46 CR and 29 ATP8 sequences. Computational analyses were performed using a number of approaches. The CR sequences were used to create a median-joining network using PhyloNet (Helgason unpublished). Bayesian phylogenies were constructed using MrBayes, both for the CR data and a dataset where ATP8 and CR sequences were combined (Huelsenbeck and Ronquist, 2001). The merged dataset consisted of the subset of samples that yielded both mtDNA sequences.

Nucleotide diversity, π; gene diversity, Hd; number of polymorphic sites, S; and Tajima’s D, Td (table 1) were calculated using DNAsp 5.0 (Librado and Rozas, 2009) on three partitions of the dataset (control region; ATP8; both markers merged). For both genetic markers, ancient DNA sequences were binned according to three main geographic areas (Europe, Siberia and Alaska plus Yukon). Given the rather limited number of sequences available from Europe, these were considered in a single statistical group, spanning a time range of 52–16,000 cal a BP. Similarly, the sequences generated for Alaskan/Yukon cave lions were grouped into a single unit which spanned a more limited time frame (22–14,000 cal a BP). Conversely, a larger set of sequences was available from Siberian specimens. We therefore further divided the Siberian data set into two time bins (62–48,000 cal a BP and 33–14,000 cal a BP) with numbers of sequences comparable to those in the other geographic bins. Pairwise FST-values were computed across Siberian and Alaskan populations to test for population differentiation, considering the two Siberian time bins separately (table 1). In order to counteract biases due to temporal asynchrony in the sampling (Depaulis et al., 2009), we complemented the FST estimation with a FT statistic (Sandoval-Castellanos, 2010), which is an FST designed to compare temporarily spaced samples. For comparing cave lions with other species of the Panthera genus, the same statistics were estimated for the extinct American lion (Panthera atrox) and the extant African lion (Panthera leo).

Table 1

Population summary statistics. Nucleotide diversity, π; gene diversity, Hd; number of polymorphic sites, S; and Tajima’s D, Td. Standard deviations of nucleotide diversity and gene diversity estimates are provided between parentheses. CR: Control Region. Time Period is given in thousands (K) of calendar years before present. *: sample EE26 was excluded for estimating the different summary statistics as the sequence recovered was only partial. NA: Not determined as the number of sequences available was not sufficient.

Gene partition Population Time Period N Length (bp) π Hd Td S

CR Europe 52–16K 8 205 0.01167 (0.00218) 0.929 (0.084) -0.5409 7
Siberia 33–14K 10 206 0.00367 (0.00157) 0.551 (0.164) -1.0345 3
Siberia 62–48K* 13 206 0.00884 (0.00190) 0.872 (0.067) -1.138 8
Alaska 22–14K 9 204 0.00191 (0.00081) 0.389 (0.164) 0.1565 1
P. atrox 33–13K 3 195 0.01368 (0.00645) 0.667 (0.314) NA 4
P. leo Modern 16 201 0.00829 (0.00160) 0.833 (0.072) -0.7398 7
ATP8 Europe 52–16K 3 143 0.01399 (0.00491) 1.000 (0.272) NA 3
Siberia 33–14K 9 143 0.01049 (0.00478) 0.583 (0.183) -1.3984 6
Siberia 62–41K 14 143 0.00200 (0.00166) 0.143 (0.119) -1.4807 2
Alaska 22–14K 3 143 0 0 NA 0
P. leo Modern 3 143 0.01399 (0.00491) 1.000 (0.272) NA 3
CR + ATP8 Europe 52–16K 3 349 0.01146 (0.00403) 1.000 (0.272) NA 6
Siberia 33–14K 9 349 0.00669 (0.00233) 0.806 (0.120) -1.3594 9
Siberia 62–41K 11 349 0.00636 (0.00173) 0.873 (0.089) -1.4912 10
Alaska 22–14K 3 348 0.00192 (0.00090) 0.667 (0.314) NA 1
P. leo Modern 3 350 0.01143 (0.00381) 1.000 (0.272) NA 6

We further investigated the demographic dynamics of cave lion populations using the Bayesian phylogenetic inference software BEAST v1.6.1 (Drummond and Rambaut, 2007) on either the control region or merged partitions. No BEAST analysis was conducted on the ATP8 gene alone, mainly due to limited variation of the generated fragment. More specifically, analyses were performed under a constant size population model and under relaxed demographic models (extended Bayesian skyline and GMRF Bayesian skyride) and the respective performance of the different models was assessed through the analysis of Bayes factors (Suchard et al., 2001).

All analyses employed a Kimura-2 parameter model (i.e. HKY with equal base frequencies) with homogeneous rates (no gamma correction) because HKY was the preferred model from the analysis in Modeltest and in order to keep in line with the coalescent simulations that only can be performed under a Kimura 2-parameter model (see below). In addition to this, a relaxed-clock model (uncorrelated lognormal distribution across branches) was assumed. Analyses where a strict-clock model was considered provided similar results (data not shown), but we retained the relaxed-model since it has been proved to contain less bias (Battistuzzi et al., 2010; Drummond et al., 2006). A total of 25 million generations was performed, where the first 7.5 million steps were discarded as burn-in. Convergence and ESS values superior to the critical threshold of 200 were checked using Tracer 1.5 (Drummond and Rambaut, 2007).

Additionally, approximate Bayesian computation coupled with serial coalescent simulations was performed on the CR data set of the Siberian and Alaskan (Beringian) samples, for a closer examination of the population dynamics. In a first stage, a model choice analysis was performed to compare eight alternative population models, assuming different demographic and migration dynamics across the Bering Strait over time (figure 2). A first model showing no demographic change over time and no geographic structure across Beringia was considered as a null alternative. The second model was similar to the null model, except that it incorporated a demographic decline of a 1- to 10-fold magnitude occurring in the population between 50–25,000 cal a BP. This simulation set-up was similar to the one described in a recent study that defined statistical guidelines for improving the ability to detect demographic bottlenecks using time-series of DNA sequence data (Mourier, 2012). Model 3 contained a similar decline as in the previous one, but with a subsequent expansion in population size between 25–12,000 cal a BP. All further five population models assumed a geographic structure between Siberian and Alaskan populations (figure 2), with (models 5, 6 and 8) or without migration (models 4 and 7) between both continents. Contrary to models 5 and 8, which assume constant gene flow between populations, migration was stopped at 50–25,000 cal a BP in model 6 and no further migration was allowed up to the age of our most recent specimen. In addition, complex population histories combining geographic structure and independent demographic bottlenecks were reconstructed in models 7 and 8. Overall, the whole set of models considered allowed us to explore which demographic and migration dynamics were better supported by the genetic information available.

Figure 2 

Population models tested using serial-coalescent simulations. Letters indicate the parameters sampled from priors distributions (other than mutation rates): T, M, I, and Ne, all sampled from uniform distributions. Dashed line indicates the timing (t) of important events: sudden change in population size (priors set to 11.2-25,000 cal a BP in the most recent event of model 3, and 25-60,000 cal a BP for the oldest event in model 3 and events in models 2, 7 and 8), stop in migration (coinciding with times of demographic change), and geographic split (60-75,000 cal a BP in models 5 to 8). Arrows indicate migration (M) set with priors ranging from 0 to 1% of migrants per generation while I is the intensity of the demographic change set with a prior between 1 to 10 times the size of the previous or posterior population change. Support values are listed below each model, for the two calibration rates used. Maximal marginal likelihood is reported with an asterisk.

In a second stage, after discarding the models that received significantly lower support, a new model choice analysis was performed for testing only single-population models where Siberian and Alaskan (Beringian) samples were considered as a single population (models 1, 2 and 3 in figure 2). Following this second-stage model choice analysis, we also performed additional analyses of key parameters in the selected model, such as effective population size, as well as timing and magnitude of the demographic changes. We employed an MCMC without likelihoods algorithm in order to perform a rejection on-the-go with an acceptance threshold that would accept 1% of the simulations under a regular procedure (a rejection performed a posteriori). A total of 5 million simulations obtained under the MCMC scheme were subjected to a further rejection procedure to obtain 5,000 final simulations (0.1%). This approach yielded a posterior sample that is equivalent to what could be obtained by running half a billion simulations, all with the goal of achieving a high accuracy in the posterior distributions of the parameters. Simulations for stages 1 and 2 were performed using the program Bayesian Serial SimCoal (Anderson et al., 2005), while the software BaySICS (Sandoval-Castellanos et al., 2014) was used for the MCMC procedure and support analyses.

The generation time was set to five years (Björklund, 2003), whilst the mutation model had a transition/transversion bias of 0.8801 and a rate sampled from a prior with normal distribution with mean of 2.57.10–7 mutation/site/yr and standard deviation equal to 1.2.10–7. These parameters were estimated from the posterior distributions of the BEAST analyses described above. Since this first mutation rate was obtained through tip calibration (Ho et al., 2007), we also performed a second set of analyses using a mutation rate derived from a BEAST analysis where a calibration point at 550 ± 25 thousand years was included for the split between the Panthera spelaea and Panthera leo groups, following Barnett et al. (2009), in order to improve our estimation and reduce effects of time-dependency (Shapiro and Ho, 2014). This second approach resulted in a slightly lower estimate for the mutation rate (1.05x10–7 mutation/site/yr with a standard deviation of 2.7.10–8). Thus we performed the entire analysis also with this updated mutation rate. However, using this second lower mutation rate did not have any effect on the model choice in the coalescent simulations (see below). Priors for population sizes were uniform (100–100,000 initially and after an optimization round they were set to 50–50,000). The summary statistics that were employed comprised the number of haplotypes, number of segregating sites, nucleotide diversity and Tajima’s D for each statistical group. Additionally, the number of pairwise differences and FST values between the statistical groups were used. Finally, posterior probabilities of all models were estimated employing simple acceptance ratios and categorical regression following Beaumont (2008).

Results

From the 29 samples, 11 yielded consensus sequences for both ATP8 and CR, one yielded a partial CR sequence (CRB) and two yielded only the ATP8 sequence (table 1). Six samples were cloned due to fragmentation and two were cloned after detection of nuclear insertions. The final consensus sequences were confirmed by two sequences from independent PCR amplifications. In cases where sequences from the two amplifications differed (samples EE3, EE16 and EE26) a third round of amplification and sequencing was conducted to obtain a majority-rule consensus sequence. The identified sequence differences consisted mainly in GC → AT transitions (83%), most likely resulting from PCR nucleotide misincorporations due to post-mortem Cytosine deamination (Hofreiter et al., 2001). All blank controls were negative, although a sequence from domestic dog (Canis lupus familiaris) was observed in one sample (EE26). However, additional amplification and sequencing of this specimen yielded a Panthera spelaea sequence.

Four previously reported CR haplotypes were identified for six samples, whereas the other five represented new haplotypes (table S1). Two of these new haplotypes were collected from southern Siberia/Baikal (EE6 and EE7: Y1 and EE16: Y2), a region from where no previous samples had been analysed. The other two (EE3: Y3, EE14: Y4) came from the Urals and North-East Siberia respectively. The sample that yielded only a partial CR sequence (EE26) did not conform to any known haplotype, and is therefore referred to as Y5 (table S1).

The median-joining network showed that the CR haplotypes can be divided into two groups, henceforth referred to as haplogroups I and II (figure 4). Although it yielded only a partial sequence, the sample EE26 (41,025 ± 970 cal a BP) did show a closest relation to one of the ‘older’ haplotypes (J), and could thus be grouped within haplogroup I (figure 4). All samples younger than this date belong to haplogroup II. Assuming a binomial sampling, this shows that the frequency of haplogroup I was below 12.3% and could even suggest that haplogroup I might have disappeared completely sometime after 41,000 cal a BP.

Figure 3 

Posterior distributions of the parameters in model 3. The posterior distributions were estimated using an approximate Bayesian computation procedure and a mutation rate resulting from a tip and split calibration with BEAST v1.6.1. Prior distributions are shown as bars and dotted lines indicate 95 % credibility intervals. a) Female effective size (Nfe) of the population following the demographic bottleneck (number of haploid gene copies). b) Magnitude of the demographic change of each event as a ratio of the larger population size to the smaller one; decline (right curve) and recovery (left curve). c) Timing of the demographic decline (right curve) and recovery (left curve).

Figure 4 

Haplotypes in network and temporal graph. a) Median-joining network analysis based on the control region sequence, including samples from Barnett et al. (2009). Letters represent haplotypes given in table 1 and colors correspond to geographic region. Small crossing bars represent number of mutations between haplotypes. Since haplotype Y5 is based on an incomplete sequence it is connected with a dashed line. b) Temporal arrangement of samples (including those from the study of Barnett et al. (2009)) and their respective haplotypes. The estimated bottleneck in Beringia is shaded and delimited by the mode values. Dotted lines show 95% credibility intervals. Extended and faded timelines have been added for samples lacking a finite date.

The phylogeny based on the CR (figure S1) did overall recover the same topology as Barnett et al. (2009). The monophyly of haplogroup II observed in the CR dataset (figure S1) was also supported by the phylogeny based on the combined ATP8 and CR datasets, where it was supported by a posterior probability of 100% (figure S2).

For the CR, we found that the nucleotide diversity (π) within Siberian cave lions in the period 33–14,000 cal a BP was lower than the diversity 62–48,000 cal a BP, even though the two periods spanned similar time ranges (table 1). Additionally, the former was found to be lower than the nucleotide diversity present in modern African lions (table 1), even though the nucleotide diversity in the cave lions might have been over-estimated due to age difference among samples (Depaulis et al., 2009). We found a quite different pattern in ATP8, with minimal diversity estimates in Siberia 62–41,000 cal a BP. However, we believe that the combination of a lower mutation rate and a reduced sequence length might have resulted in poorer diversity estimates using this genetic marker. All subsequent analyses were accordingly performed using the CR dataset. Pairwise FST values of 0.07 and 0.68 were obtained between the Alaskan and Siberian populations for the two time bins 62–48,000 and 33–14,000 cal a BP, respectively. The pairwise FST value among the Siberian time bins was between 0.3 to 0.6 but the FT values dropped to 0.04–0.15 (p<0.001), depending on the Ne. Significant results were also obtained between old Siberian samples and Alaskan/Yukon samples (FT = 0.01–0.13; p = 0.01–0.047). The FT being positive and significant rejects a null hypothesis that the changes were due to genetic drift in a population that is panmictic and demographically constant through time. This could be due to demographic changes, isolation between subpopulations, or adaptation.

We evaluated different demographic models implemented in BEAST (constant size, Bayesian skyline and GMRF Bayesian skyride) and tested for their difference in likelihood using Bayes factors (table 2). Regardless of the dataset considered, neither skyline nor skyride models showed better support than constant size models. It should be noted that the Bayes factors rely on harmonic-mean estimators of the marginal likelihood, which is an approach that may be biased towards parameter-rich models (Baele et al., 2012; Baele et al., 2013). However, the congruence among datasets as well as our assays by ABC seems to indicate that the lack of resolution among models is due to limited sequence length and sampling rather than the method. One possible way to overcome this would be to include additional loci (Ho and Shapiro, 2011).

Table 2

Support for different population demographic models using BEAST and tip calibration. The marginal log likelihood of each model, as well as Bayes Factors (logBF) between alternative models and best models are shown. Different sets of sequences originating from restricted or large geographic regions were considered in order to test for population structure as a possible confounding factor.

Populations Model log Likelihood logBF

Siberia Constant -375.5 (0.144) -0.355
Siberia Skyline -376.0 (0.142) -0.553
Siberia Skyride -374.7 (0.121) *
Siberia + Alaska Constant -385.3 (0.145) -0.258
Siberia + Alaska Skyline -384.7 (0.133) *
Siberia + Alaska Skyride -385.0 (0.135) -0.099
Siberia + Alaska + Europe Constant -434.4 (0.183) -0.433
Siberia + Alaska + Europe Skyline -433.4 (0.177) *
Siberia + Alaska + Europe Skyride -433.7 (0.189) -0.128

We also performed extensive serial coalescent simulations of eight different population models, ranging from the simple null scenario of constant size in a single deme to different situations with multi-demes experiencing some demographic shift, and in some cases also connected through migration (figure 2). We used categorical regression within an approximate Bayesian computation framework to estimate the posterior probability of all models. Only the model where Siberian and Alaskan lions were grouped in a single deme and experienced an instantaneous demographic decline followed by a recovery (model 3) received consistently more support than the null scenario (figure 2) with a Bayes factor of 2.0–4.0 in a direct comparison; this held true regardless of the mutation rate priors considered. The same model without a recovery (model 2) had less support, similar to the null scenario (model 1). Importantly, none of the other more complex population models investigated received even marginal support.

The second set of simulations was performed with the Siberian and Alaskan (Beringian) samples as a single population with independent population sizes. In the posterior sample of those simulations, 86–91% of them supported a demographic decline in the time interval 50–25,000 cal a BP, and 75–80% of the simulations gave support for the same decline followed by a recovery in the time interval 25–12,000 cal a BP. This suggests that the aDNA sequences recovered in Siberia and Alaska better fit a model of a single population that experienced a demographic bottleneck. Based on the posterior distributions of the final analyses (figure 3), the bottleneck likely started around 47,000 cal a BP (mode = 46,500; 95% credibility interval (C.I. = 28,482–54,754) and was followed by a recovery around 18,000 cal a BP (mode = 18,000; 95% C.I. = 12,252–24,507) to a post-bottleneck effective population size of approximately 1,000 individuals (figure 3 & table 3).

Table 3

Parameter estimations of model 3 using ABC. Estimations of: effective population size after the bottleneck (Ne1), magnitude of expansion (Ne1/Ne2), magnitude of decline (Ne3/Ne2), time (cal a BP) to the end of the decline (T1), and time to the beginning of the decline (T2).

Mode Median 95% C.I.

Ne1 1000 1476 129.64 8658.4
Ne1/Ne2 0.15 0.2796 0.0432 0.9933
Ne3/Ne2 0.05 0.0633 0.0147 0.4392
T1 18000 17481 12252 24507
T2 46500 41603 28482 54754

Discussion

The use of ancient DNA has in many ways revolutionised our understanding of past population events in many species, and with new technologies the ways to retrieve DNA from fossil remains are currently improving. But there still exist several limitations, mainly depending on available remains and their preservation. With these restrictions in mind the amount of accumulated sequence data is still modest in general and any additions are valuable, especially from extinct species. Bearing in mind that our results are based on short sequences from a single locus (the mitochondrial genome), and thus should be interpreted with caution, the additional samples analysed here nonetheless allowed us to obtain an improved resolution of the population dynamics in the cave lion during the Late Pleistocene. Moreover, the use of coalescent simulations analysed in an ABC framework enabled a formal evaluation of the hypothesis of the population decline in Beringia proposed by Barnett et al. (2009). Future analyses including genomic data will likely provide a more powerful approach to further address this hypothesis.

As shown by the temporal distribution of control region haplotypes, a marked decline in genetic diversity indeed seems to characterize the population dynamics of Late Pleistocene cave lions in Beringia (figure 4). The approximate Bayesian computation analyses supported a model where the Beringian population experienced a bottleneck event starting approximately 47,000 cal a BP (95% C.I. = 28,482–54,754) which was then followed by a recovery (model 3). This bottleneck was probably also the cause behind the eventual demise of haplogroup I (figures 3 and 4). Although the effects of the bottleneck are seen most clearly in Beringia, it should be noted that all haplotypes observed in samples from other regions that are younger than 41,000 cal a BP also belong exclusively to haplogroup II (figure 4), implying that the genetic decline might have affected these areas as well. The star-like pattern observed in haplogroup II indicates that the bottleneck was followed by a re-expansion in population size, which is also consistent with the large geographical distribution of cave lion fossils dated to between ca. 24–13,000 cal a BP (Stuart and Lister, 2011). According to the proposed bottleneck model, the timing of the re-expansion is quite consistent with this period in time (figure 3).

The timing of the bottleneck estimated from the genetic data also overlaps with a pronounced gap in the fossil record of cave lions in Siberia at ca. 40–35 cal a BP (Stuart and Lister, 2011). Whether this merely implies a reduced population size or a total but temporary extinction in most parts of the range cannot be addressed without analysis of further samples, or more extensive ancient DNA data from the existing samples (e.g. nuclear SNPs (Mourier, 2012)).

In Western Europe, there is a second gap in the fossil record between approximately 30,000 and 19,000 cal a BP (Stuart and Lister, 2011), which corresponds roughly in time with the Last Glacial Maximum (LGM). It has therefore been proposed that cave lions disappeared from Western Europe during the LGM (Stuart and Lister, 2011). Unfortunately, the small sample size from Western Europe in this study precludes any genetic assessment of this hypothesis.

In the Beringian cave lions the start of the bottleneck is approximately consistent with the timing proposed by Barnett et al. (2009) and falls within marine isotope stage (MIS) 3 (ca. 60–30,000 cal a BP), which was a period characterised by repeated climatic fluctuations on a millennial and multi-millennial scale (Grootes et al., 2001; Voelker, 2002). For a large carnivore such as the cave lion, it is possible that climate-induced changes could have had a strong impact, for instance on the availability of prey.

There are several other large mammals in the Beringia region that also display patterns of demographic decline during MIS 3, many of which constitute potential prey species of the cave lion. Mammoths (Mammuthus primigenius) experienced a genetic decline in Beringia somewhere around 44,000 cal a BP, as an entire mitochondrial DNA clade vanished around this time (Barnes et al., 2007; Gilbert et al., 2008). In steppe bison, a severe demographic reduction started some 44–34,000 cal a BP (Drummond et al., 2005; Shapiro et al., 2004). Ancient DNA analysis combined with radiocarbon dating on brown bears (Ursus arctos) has identified an extinction of a genetically unique population in eastern Beringia, which was followed by an extended gap in the fossil record between 40–25,000 cal a BP (Barnes et al., 2002). A similar gap is seen for Beringian horses (Equus sp.) at about 40–33,000 cal a BP, including the extinction of a small-bodied hemionid species at about 35,000 cal a BP (Guthrie, 2003). Finally, a study on muskox (Ovibos moschatus) has revealed a loss of a genetically distinct clade in North-East Siberia sometime after 47,000 cal a BP (Campos et al., 2010b).

Although not entirely overlapping, the numerous declines and losses of genetic variation, as well as the differing ecologies of the affected species, implies that major changes in the environment of Beringia may have occurred during this time. A potential change could have been rapid climate shifts that would have had a considerable impact on the steppe-tundra habitat in Beringia, and which also would agree with the pronounced temperature variation recorded throughout MIS 3 (Grootes et al., 2001; Voelker, 2002). One specific climate shift suggested by Barnett et al. (2009) as a possible cause of the bottleneck in cave lions was Heinrich event 5a, a cold stadial occurring around 55–52,000 cal a BP (Heinrich, 1988; Hemming, 2004). Although this event cannot be ruled out as an underlying cause, the wide time range we recovered for the estimated start of the bottleneck lends no certainty to this explanation. Moreover, our analyses suggested a prolonged bottleneck that lasted several tens of thousands of years. This could indicate that the bottleneck was caused by a long-term change in the cave lion’s environment rather than a short-term fluctuation in climate, such as a Heinrich event.

The direct causes behind the subsequent re-expansion of cave lions in Beringia remain elusive as well, but might have followed increasing prey abundance. Recent genetic analyses suggest demographic expansions in Eurasian horse and muskox approximately 34,000 cal a BP, which could indicate a subsequent recovery in these two species at the end of MIS 3, possibly in conjunction with an increase in habitat availability (Lorenzen et al., 2011). Both the paleontological record and genetic data suggest a demographic recovery in the cave lion sometime around 25–14,000 cal a BP. Although the demographic recovery in the cave lion thus appears to have occurred later than the recoveries in horse and muskox, it seems likely that more data is required to exclude that these events were somehow connected.

Interestingly, several of the species that suffered demographic declines along with cave lions during MIS 3 also went through severe declines or even extinctions in Beringia at the very end of the Pleistocene. Mammoths and muskox were confined to small arctic refugia (Campos et al., 2010b; Lister and Stuart, 2008), while horses and steppe bison disappeared entirely (Guthrie, 2003; Shapiro et al., 2004; Stuart and Lister, 2011). To what extent the declines during MIS 3 contributed to the final extinction process in these species has not yet been fully explored. However, a lowered genetic diversity might have rendered these species more vulnerable to the ecological changes that affected Beringia at the Pleistocene-Holocene boundary.

Another potential cause for population decline in carnivores is inter-specific competition. For example, in the case of the brown bear, the local extinction in eastern Beringia around 40,000 cal a BP was more or less synchronous with the arrival of the short-faced bear (Arctodus simus) from southern North America (Barnes et al., 2002). The two bear species did co-exist for a long time period prior to this, but the hypercarnivorous short-faced bear may have negatively impacted Beringian cave lion as well as scimitar-tooth cat (Homotherium serum) populations. Despite an extensive collection of radiocarbon dates from Alaska and the Yukon, there are very few cave lions recognised from MIS 3 (Stuart and Lister, 2011), a time period that also saw the probable extinction of H. serum (latest date of ca. 41,000 cal a BP) (Fox-Dobbs et al., 2008). The claim that Homotherium was present in the European Late Pleistocene, (Reumer et al., 2003) thereby overlapping with P. spelaea, needs confirmation, but both taxa occurred at several Middle Pleistocene European sites, indicating that niche partitioning was possible (Anton et al., 2005). However, the Old World species H. latidens may not have had the same dietary and habitat requirements as the new world H. serum. Stable isotope analysis shows that both scimitar-tooth cats and cave lions in Beringia appear to mainly have been either generalists or horse/bison specialists (Fox-Dobbs et al., 2008). Local extinction of scimitar-tooth cat in eastern Beringia, whether due to the short-faced bear or underlying habitat changes, may therefore have facilitated an expansion of cave lions from their western range during MIS3/2 and subsequent bottleneck recovery.

Extinctions and population declines in large mammals can often be associated with human encroachment and overexploitation, and this has also been suggested as a contributing factor to the end-Pleistocene extinction of the Eurasian and American megafauna (e.g. (Barnosky et al., 2004). The inferred start of the demographic decline in the cave lion at ca 47,000 cal a BP (95% C.I. = 28,482–54,754) overlaps with the first appearance of anatomically modern humans in parts of the cave lion’s range, such as Europe and southern Siberia (Fu et al., 2014; Goebel, 1999; Higham et al., 2011). In Beringia, the earliest well-supported evidence of human settlement in north-eastern Siberia is dated to ca. 32,000 cal a BP (Pitulko et al., 2004). Thus, although expansions of modern humans may have contributed to the decline in some parts of the cave lion’s distribution, such as Europe and southern Siberia, we find it less likely that the onset of the demographic bottleneck in Beringian cave lions was the result of anthropogenic factors.

Conclusions

In conclusion, this study indicates a demographic decline in the cave lion starting around 47,000 cal a BP, which was then followed by a recovery at some point roughly 18,000 cal a BP. Moreover, it seems that several other taxa in Beringia, including some of the cave lion’s potential prey species, may have gone through similar declines. The mechanism behind this pattern remains uncertain, although the declines all seem to take place during MIS 3 and in Beringia some pre-date the arrival of anatomically modern humans. Further research is required to investigate other potential causes of the declines, such as climate-driven habitat change and interaction among species in this region. A combined approach, considering both genetics and paleoecology, should allow a better understanding of past population dynamics.

Acknowledgements

We are grateful to Marina Sotnikova, Doris Nagel, Gennady Baryshnikov, Aleksandr Matuzko, Tatiana Kuznetsova, Alexander A. Shchetnikov, Dustin White and Sergey Vartanyan for providing samples for the DNA analysis. This study was financed through support from the Swedish Research Council and Formas through the FP6 BiodivERsA ERA-NET program CLIMIGRATE. We also thank Agnar Helgason for making the software PhyloNet available, as well as Pia Eldenäs, Martin Irestedt, Bodil Cronholm and Keyvan Mirbakhsh for assistance in the laboratory. Ludovic Orlando was supported by a Marie-Curie Career Integration Grant CIG-293845. Edson Sandoval Castellanos acknowledges support from the strategic research programme EkoKlim. Adrian Lister and Anthony Stuart gratefully acknowledge the support of the Natural Environment Research Council (Grant NE/G00188X/1), and R.B. acknowledges the support of a Marie Curie Intra-European Fellowship (FP7-PEOPLE-2011-IEF-298820).

Competing interests

The authors declare no competing financial, non-financial or other interests in relation to this manuscript.

Availability of supporting data

Sequence data produced in this study have been submitted to GenBank (accession numbers JX486724-JX486734 and KP280033-KP280046). Data from previously published studies are available from GenBank (accession numbers DQ899900-DQ899945).