Background

The extinct cave lion (Panthera spelaea) was an integral component of the late Pleistocene Holarctic ecosystem, occupying the position of apex predator (, , , ) alongside the scimitar cat (Homotherium sp.) (). The cave lion is known from plentiful remains preserved in the karstic cave systems of Europe and the permafrost sediments of Beringia (present day Siberia, Alaska and the Yukon). Despite this, P. spelaea has had an interesting history of taxonomic revision, having been variously considered a highly divergent population of the modern lion (), a subspecies of the modern lion (), and a full species in its own right (). There is currently no consensus on the taxonomic status of the cave lion.

Cave lions were significantly larger than their modern lion counterparts, and exhibit a unique cranial morphology that has a mosaic of characters present in the lion and tiger (, , , ). Evidence from Pleistocene art demonstrates that the cave lion differed in external morphology from modern lions (). Male cave lions did not possess a mane, a notable secondary sexual character in the modern species, and this is likely to have had considerable effect on the ethology of the species (, ). Previous work has also demonstrated that the cave lion sensu stricto had an enormous range (), from western Europe to eastern Beringia, and were coeval with lion populations in southern North America () and Africa (Panthera leo; ). Despite dominating the Holarctic mammoth steppe for most of the late Pleistocene, the cave lion went extinct nearly simultaneously across its range, with terminal dates in Europe, Siberia, and Alaska all close to 14,000 cal BP. Panthera atrox disappeared only slightly later in southern North America ().

Ancient DNA (aDNA) studies based on partial fragments of the mitochondrial control region, ATP8, and cytochrome b genes from P. spelaea have not fully resolved the degree of separation between the various lineages of lion-like cats (, , ). In particular, although prior analyses recovered spelaea and leo as sister taxa, these studies have relied on the Middle Pleistocene appearance of the ancestral cave lion (Panthera [leo] fossilis) to calibrate the age of the split for molecular dating analyses, without estimating the timing of the split directly. Here, we present the complete mitochondrial genomes of two cave lion specimens and, using a Bayesian phylogenetic approach, we analyse these sequences in combination with a strictly vetted set of published mitochondrial genomes. We use these to infer the evolutionary timescale of lion-like cats, calibrated using recently discovered pantherine fossils, and discuss the species status of P. spelaea.

Materials and Methods

Sample Discovery

Sample YG 401.410, a humerus (Figure 1, Table 1), was recovered from the Quartz Creek site in Yukon, Canada (Figure 2), on 24 July 2010. It has subsequently been stored in the Government of Yukon’s Palaeontology Program collection in Whitehorse. The specimen was subsampled in June 2013 by GZ and sent to the Centre for GeoGenetics (University of Copenhagen, Denmark) for processing. Hair sample F-2678/70 was found at the right bank of Malyi Anyui river (Chukotka, Russia) during the summer of 2008 (). A small bundle of hair was subsampled by IVK and sent to the UCSC Paleogenomics Lab (Santa Cruz, USA) for processing.

Figure 1 

A) Four views of a humerus of cave lion (Panthera spelaea), specimen YG 401.410. This specimen was recovered from Quartz Creek in the Yukon Territory, Canada. The scale bar is 5 cm. B) Four views of cave lion (Panthera spelaea) hair bolus, sample F-2678/70, contains guard hair but is mainly represented by thick underfur of tightly packed, wavy fur hairs. This specimen was recovered from Malyi Anyui river in Chukotka, Russia.

Table 1

Metric data from humerus specimen YG 401.410. Measurements to nearest 0.1 mm taken using digital calipers.

SpecimenSideMinimum breadth of diaphysis (mm)Maximum depth of diaphysis (mm)Maximum breadth of distal end (mm)Minimum anteroposterior diameter of articulating surface for ulna (mm)Notes
YG 401.410Left28.652.186.129.1Missing proximal end from pectoral ridge
Figure 2 

A map of the approximate distributions of late Pleistocene lions and the provenance of samples used in this study. Red indicates the maximal range of Panthera spelaea; blue the maximal range of Panthera atrox; and green the maximal range of Panthera leo leo/Panthera leo persica. Stars show approximate locations of lion samples. The insets show details of the modern boundaries of Yukon Territory, Canada, and Chukotka, Russia, along with regional settlements. Sample YG 401.410 was found at Quartz Creek (Yukon Territory, Canada), 63.49N, 139.27W. Sample F-2678/70 was found at the right bank of Malyi Anyui River (Chukotka, Russia), 68.18N, 161.44E.

Radiocarbon Dating

A section of bone from sample YG 401.410 was delivered to Stafford Research LLC (University of Aarhus, Denmark). Samples of crushed bone were decalcified and washed, treated with 0.05 N NaOH overnight to remove humics, soaked in 0.1 N HCl, gelatinized at 60°C at pH 2, and ultrafiltered at 30 kDa. The purified collagen was then graphitised and analysed at the W.M. Keck Carbon Cycle Accelerator Mass Spectrometry (AMS) Laboratory, University of California Irvine, according to standard protocol (, , ).

Extraction and DNA Amplification

Sample YG 401.410

Samples of cortical bone were taken (approx. 1 cm3) using a Dremel powertool and reduced to powder in a Mikrodismembrator. DNA extraction was performed as described by Orlando et al. () in a dedicated ancient DNA laboratory at the Centre for GeoGenetics in Copenhagen, in parallel with negative extraction controls.

The DNA extract and negative control were then built into genomic libraries using the NEB E6070 kit (New England Biolabs), following a protocol slightly modified from that by Vilstrup et al. (). Briefly, extract (30 µL) was end-repaired and then passed through a MinElute column (Qiagen). The collected flow-through was then adapter-ligated and passed through a QiaQuick column (Qiagen). Adapter fill-in reaction was then performed on the flowthrough, before final incubation at 37°C (30 min) followed by inactivation overnight at –20°C.

We then amplified the DNA in a 50 µL reaction, using 25 µL of library for 12 cycles under the following reaction conditions. Final concentrations were 1.25 U AccuPrime™ Pfx DNA Polymerase (Invitrogen), 1× AccuPrime™ Pfx reaction mix (Invitrogen), 0.4 mg/mL BSA, 120 nM primer in PE, and 120 nM of a multiplexing indexing primer containing a unique 6-nucleotide index code (Illumina).

PCR cycling conditions consisted of an initial denaturation step at 95°C for 2 min, followed by 12 cycles of 95°C denaturation for 15 s, 60°C annealing for 30 s, and 68°C extension for 30 s. A final extension step at 68°C for 7 min was also included. Amplified libraries were first checked for presence of DNA on a 2% agarose gel before purification using the QIAquick column system (Qiagen) and quantification on an Agilent 2100 BioAnalyzer. Quantified libraries were communally pooled in equimolar ratios and sequenced as single-end reads (100 bp) on an Illumina HiSeq2000 platform at the Danish National High-Throughput Sequencing Centre.

Sample F-2678/70

DNA extraction from hair sample F-2678/70 was performed in a dedicated, sterile, facility at UC Santa Cruz, using standard protocols for ancient DNA (). The extraction followed the protocol described by Dabney et al (), with the modifications suggested by Campos and Gilbert (). The Illumina DNA sequencing library was built following Meyer and Kircher (). Between each step, the libraries were cleaned using Sera-Mag SPRI SpeedBeads (ThermoScientific) in 18% PEG-8000. The resulting libraries were sequenced in a MiSeq Illumina sequencer with v3 kits at UCSC Paleogenomics Lab.

Sequence Processing

Sample YG401.410

The dataset consisted of 3,444,248 single-end sequences that were cleaned of adapter and low-quality sequences using AdapterRemoval v1.2-GG1 (). In order to obtain the mitochondrial sequence, we mapped the resulting 3,429,186 cleaned reads against a published mitochondrial genome from Panthera leo persica (GenBank accession JQ904290.1) using bwa (). In order to take into account the circularity of the mitochondrial genome and to recover the reads mapping at the edges of the extremes of the reference sequence, we added the first 100 bases to the end. Clonality was removed with picard-tools v1.92 (). The mapping was analysed for damage at the extremes with MapDamage v2.0 () and the quality of the damaged bases was rescaled with MapDamage. Realignment was performed with GenomeAnalysisTK () on the rescaled mapping file in order to call SNPs with samtools v0.1.18 and bcftools () with a minimum coverage of 8 and minimum quality of 30.

Nuclear copies of mitochondrial genes (numts) are known to occur in multiple members of the cat family Felidae. Some phylogenetic studies of felids have been compromised by the inclusion of numt sequences in alignment (). Numts were anticipated to be a particular problem in genomes reassembled from high-throughput sequencing technologies, which involve short fragment lengths and are unable to preferentially target cytoplasmic copies (e.g. by PCR primer design).

In order to obtain a robust, non-chimaeric consensus sequence with GenomeAnalysisTK, we only used those SNPs supported by at least 2/3 of the reads mapping to that position. We used the resulting consensus sequence as a reference for a second mapping round, using the same mapping and consensus strategy as before. This was repeated three more times to make a total of four mapping rounds. The final consensus sequence has an average depth of 9.53× covering 89.5% of the mitogenome.

F-2678/70

The dataset consisted of 6,683,556 paired-end reads. Initial processing of these reads consisted of using AdapterRemoval v1.2-GG1 () for the cleaning step, followed by merging reads including a minimum overlap of 10 base-pairs between forward and reverse reads using SeqPrep (http://github.com/jstjohn/SeqPrep). A total of 6,123,696 merged reads were obtained. These merged reads were mapped against the Panthera leo persica mitochondrial genome (GenBank accession NC_018053.1) using MIA (http://github.com/udo-stenzel/mapping-iterative-assembler), a reference based, iterative, short-fragment assembler that accepts circular genomes as reference sequences. A total of 9,582 reads mapped to the reference mitochondrial genome. The consensus sequence was called from the resulting output file, with each base having a minimum depth of coverage of 3× and 2/3 base agreement. The final assembly had an average depth of 28.16× covering 99.62% of the mitochondrial genome.

Phylogenetic Analysis

Complete mitochondrial genomes of lion (Panthera leo), leopard (Panthera pardus), jaguar (Panthera onca), snow leopard (Panthera uncia), tiger (Panthera tigris), clouded leopard (Neofelis nebulosa), and domestic cat (Felis sylvestris catus) were downloaded from Genbank (Table 2). Sequences of all 37 mitochondrial genes were extracted from these genomes.

Table 2

Mitochondrial genomes of pantherine cats analysed in this study.

TaxonCommon nameGenbank accessionReference
Panthera leo persicaAsian lionNC018053()
Panthera leo leoAfrican lionKF776494()
Panthera leoLionKP202262Unpublished, GenBank
Panthera pardusLeopardNC010641()
Panthera pardus japonensisNorth Chinese leopardKJ866876()
Panthera oncaJaguarNC022842Unpublished, Genbank
Panthera tigris sumatraeSumatran tigerJF357970()
Panthera tigris tigrisBengal tigerJF357968()
Panthera unciaSnow leopardNC010638(, )
Neofelis nebulosaClouded leopardNC008450()
Felis sylvestris catusDomestic catFCU20753()
Felis sylvestris catus NUMTNuclear pseudogeneFCU20754(, )
Panthera tigris NUMTNuclear pseudogeneDQ151551()
Panthera spelaea YG401.410Cave lionKX258451This study
Panthera spelaea F-2678/70Cave lionKX258452This study

We estimated the evolutionary relationships among pantherine cats using maximum likelihood in RAxML v8.0.14 (). The data set was partitioned into five subsets: the three codon positions of the 13 protein-coding genes; the 2 rRNA genes, and the 22 tRNA genes. A separate GTR+G model of nucleotide substitution was assigned to each subset of the data. The maximum-likelihood tree was estimated using 10 random starts. Node support was estimated using 1000 bootstrap replicates.

To co-estimate the phylogenetic relationships and evolutionary timescale, we analysed the mitochondrial genome sequences using a Bayesian phylogenetic approach in BEAST 1.8.2 (). The Bayesian information criterion was used to select the best-fitting model of nucleotide substitution for each data subset.

We compared two tree priors (Yule and birth-death) using Bayes factors, based on marginal likelihoods calculated using the stepping-stone estimator (). To account for the potential presence of rate variation across lineages, we also used Bayes factors to compare the strict clock against the uncorrelated lognormal relaxed clock (). We included relative-rate parameters to allow each subset of the data to have a different evolutionary rate.

To calibrate our phylogenetic estimates of divergence times, we included age constraints based on the fossil record. Our calibrations were based on the stem snow leopard Panthera blytheae (), stem tiger Panthera zdanskyi (), and stem pantherine Panthera paleosinensis (). These fossils were used to specify uniform priors on the ages of corresponding nodes in the tree (). Previous estimates of the evolutionary timescale of cave lions were calibrated using fossil evidence from P. fossilis (, ), but the exact placement of this taxon is unclear. For example, incorrect assignment of a stem taxon to a crown clade can lead to overestimation of divergence times. Alternatively, divergence times can be underestimated when a taxon belonging to the crown group is erroneously assigned to the stem lineage. By using other fossil calibrations, we were able to test whether the split between cave lion and modern lion coincided with the existence of P. fossilis.

Posterior distributions of parameters, including the tree and divergence times, were estimated using Markov chain Monte Carlo sampling. Samples were drawn every 2000 steps over a total of 20,000,000 steps. Convergence was checked using two independent chains, and the resulting samples were combined. Sufficient sampling was confirmed by inspection of effective sample sizes of parameters, which were all greater than 200.

Results

Radiocarbon dating

Sample YG401.410

We obtained two separate dating results for sample YG401.410. The dates were within three hundred 14C years of each other, with overlapping 95% confidence intervals. Collagen yield was substantial (Table 3), with δ13C and δ15N values appropriate for the trophic level of the species (). Interestingly, the radiocarbon date for this specimen falls within a noticeable chronological gap for Eastern Beringian lions, suggesting that some of the absences observed by Stuart and Lister () may be resolved with further sampling.

Table 3

Results of two AMS radiocarbon analyses of sample YG401.410 with associated stable isotope and chemical analysis values.

UCIAMSUCIAMS-142833UCIAMS-143525
Fraction Modern0.0243 ± 0.00060.0234 ± 0.0006
D14C (‰)–975.7 ± 0.6–976.6 ± 0.6
14C Age BP29,860 ± 21030,160 ± 220
>30 kD Collagen yield(%)8.1N/A
δ15N (‰)8N/A
δ13C (‰)–18.2N/A
%N16.3N/A
%C44.7N/A
C/N (wt%/wt%)2.75N/A
C/N (atomic)3.21N/A

Sample F-2678/70

The hair sample has previously been radiocarbon dated by Kirillova et al. (). The reported uncalibrated AMS date of 28,690 ± 130 was much younger than dates associated with bone and claw from the same individual, which could have resulted from incomplete removal of modern carbon. For a full discussion see Kirillova et al. ().

Phylogeny and divergence times

Our initial estimates of the phylogeny gave unexpected placements for some of the taxa (data not shown). Inclusion of the mitochondrial genomes from the African lion () and Asian lion () produced a tree with P. spelaea as the sister taxon of P. l. persica, in contradiction with previous studies (, , ). Analysis of the ATP8 gene, which has been characterised for both cytoplasmic and nuclear origin () in Panthera cats, revealed the presence of a numt sequence within the mitochondrial genome published by Ma and Wang (). Given that this mitochondrial genome sequence was assembled from published nuclear genome data (), it is likely to include a significant proportion of nuclear-derived sequence. Therefore, this mitochondrial genome was excluded from further analyses.

The evolutionary relationships estimated using maximum-likelihood and Bayesian methods were congruent (Figure 3), with all nodes being strongly supported. The two samples of P. spelaea group together as the sister lineage to P. leo. Our estimates of divergence times largely overlap with previous estimates (, ). However, the estimated split between P. leo and P. spelaea at 1.89 Ma (95% CI: 1.23–2.93 Ma) is considerably older than the first appearance of Panthera fossilis, which has been used as a calibrating node in previous studies (, ).

Figure 3 

Bayesian estimate of the pantherine phylogeny based on mitochondrial genomes. The tree is drawn to a timescale, given in millions of years. Bars at nodes give 95% credibility intervals of node-age estimates. Values at nodes denote posterior probabilities and likelihood bootstrap support. The position of the root was inferred by using the domestic cat (Felis sylvestris catus) as the outgroup.

Discussion

Evolutionary history of the lion-like cats

Despite their global range and continued dominance of ecosystems in Africa, and until recently in Asia, the lion-like cats have left a confusing fossil trail. Remains of pantherine felids have been found in fossil beds dating to 3.46 million years ago at Laetoli in Tanzania (), with recognisably leonine fossils known from Olduvai II at 1.4–1.2 million years ago (). Lion fossils only become relatively abundant during the Middle Pleistocene, with the appearance of Panthera fossilis. This taxon is considered an ancestral form of Panthera spelaea and provides a minimum age for the separation between the spelaea and leo lineages. P. fossilis is known from MIS 17-15 (680-600 ka) from European sites such as Mosbach (Germany), Pakefield (UK) and Isernia (Italy) (, , , ). Given our divergence-time estimates (Figure 3), it would appear that P. fossilis must be already on the branch leading towards P. spelaea rather than close to the split. In view of this finding, we recommend that P. fossilis should only be used to provide a minimum age constraint for the split between the leo and spelaea lineages.

Dating the divergence between spelaea and leo

Much of the discussion of the taxonomic position of the cave lion has revolved around its degree of separation from modern lion populations. Although some authors have aligned the spelaea and atrox lineages with the tiger (, ) and with the jaguar (, ), most have realised its close connection to the modern lion (, , , , , ).

Previous genetic studies (, ) have used the first appearance of the ancestral cave lion (Panthera fossilis) (, ) to calibrate estimates of the pantherine evolutionary timescale. This study represents the first attempt to identify the divergence bounds between spelaea and leo without recourse to P. fossilis as a calibration. The identification of this split within the Early Pleistocene at 1.89 Ma, rather than the Middle Pleistocene, shows that the modern lion and cave lion lineages represent substantially distinct taxa. A caveat to this is that the estimates rely strongly on the fossil calibrations used. If reanalysis later shows that these fossils have been ascribed to the wrong lineages or are re-dated to different periods then their utility in the analysis will be compromised. Comparison with other recent pantherines (Figure 3) demonstrates that the degree of mitochondrial divergence is considerably greater than that found between well-defined subspecies in modern lion (P. leo) (), leopard (P. pardus) (), or tiger (P. tigris) (). The estimated divergence time between P. leo and P. spelaea is also greater than that between the two newly recognised species of clouded leopard, Neofelis nebulosa and N. diardi, which has been estimated at 1.41 Ma ().

Mitochondrial data and the considerable morphological differences () support the recognition of the cave lion as a full species: Panthera spelaea. Data from the nuclear genome will allow further testing of this proposal.

Conclusions

Our analyses of mitochondrial genome sequences reveal that the Middle Pleistocene Panthera fossilis is likely to represent a form already on the spelaea lineage. Its appearance in the fossil record demonstrates the initial spread of the ancestral cave lion form into Eurasia. Furthermore, our study has provided an estimate of 1.89 Ma (95% CI: 1.23–2.93 Ma) for the split between the lineages leading to the cave lion and modern lion. This molecular estimate appreciably antedates the appearance of P. fossilis, and provides further evidence that the cave lion was distinct enough to be considered a species in its own right.

Availability of Supporting Data

Sequence data produced for this study have been uploaded to GenBank, with accession numbers KX258451 and KX258452. Raw data have been uploaded to the Sequence Read Archive at NCBI under study accession number SRP075782.