Mitogenomics of the Extinct Cave Lion, Panthera spelaea (Goldfuss, 1810), Resolve its Position within the Panthera Cats

* Centre for GeoGenetics, Natural History Museum of Denmark, University of Copenhagen, Øster Voldgade 5-7, Copenhagen, Denmark † Department of Ecology and Evolutionary Biology, University of California Santa Cruz, Santa Cruz, CA 95064, USA. ‡ School of Life and Environmental Sciences, University of Sydney, Sydney NSW 2006, Australia § Department of Tourism and Culture, Government of Yukon, PO Box 2703, Whitehorse, Yukon, Canada ǁ Department of Biological and Environmental Sciences, Qatar University, Doha, Qatar ¶ Ice Age Museum, National Alliance of Shidlovskiy “Ice Age”, 119 bld, Mira pr, Moscow, 129223, Russia ** Palaeogenomics and Bio-Archaeology Research Network, Research Laboratory for Archaeology, Dyson Perrins Building, South Parks Road, Oxford, OX1 3QY, UK †† Trace and Environmental DNA Laboratory, Dept of Environment and Agriculture, Curtin University, Perth, WA 6102, Australia Corresponding author: Ross Barnett (drrossbarnett@gmail.com) RESEARCH PAPER


Background
The extinct cave lion (Panthera spelaea) was an integral component of the late Pleistocene Holarctic ecosystem, occupying the position of apex predator (Barnett et al., 2009, Antón et al., 2005, Yamaguchi et al., 2004, Bocherens et al., 2011 alongside the scimitar cat (Homotherium sp.) (Barnett, 2014). 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 (Dawkins et al., 1866), a subspecies of the modern lion (Kurtén, 1985), and a full species in its own right (Sotnikova and Nikolskiy, 2006). 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 (Sotnikova and Nikolskiy, 2006, Hemmer, 2011, Vereshchagin, 1971, Gromova, 1932. Evidence from Pleistocene art demonstrates that the cave lion differed in external morphology from modern lions (Packer and Clotte, 2000). 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 (Yamaguchi et al., 2004, Nagel et al., 2003. Previous work has also demonstrated that the cave lion sensu stricto had an enormous range (Barnett et al., 2009), from western Europe to eastern Beringia, and were coeval with lion populations in southern North America (Panthera atrox; Montellano-Ballesteros and Carbot-Chanona, 2009) and Africa (Panthera leo; Bougariane et al., 2010). 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 (Stuart and Lister, 2010). 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 (Burger et al., 2004, Barnett et al., 2009, Ersmark et al., 2015. 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.

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 (Kirillova et al., 2015). A small bundle of hair was subsampled by IVK and sent to the UCSC Paleogenomics Lab (Santa Cruz, USA) for processing.

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 (Stafford et al., 1988, Waters et al., 2015, Beaumont et al., 2010. Sample YG 401.410 Samples of cortical bone were taken (approx. 1 cm 3 ) using a Dremel powertool and reduced to powder in a Mikrodismembrator. DNA extraction was performed as described by Orlando et al. (2013) in a dedicated ancient DNA laboratory at the Centre for GeoGenetics in Copenhagen, in parallel with negative extraction controls.

Extraction and DNA Amplification
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. (2013). 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 (Cooper and Poinar, 2000). The extraction followed the protocol described by Dabney et al (2013), with the modifications suggested by Campos and Gilbert (2012). The Illumina DNA sequencing library was built following Meyer and Kircher (2010). Between each step, the libraries were cleaned using Sera-Mag SPRI SpeedBeads (ThermoScientific) in  The dataset consisted of 3,444,248 single-end sequences that were cleaned of adapter and low-quality sequences using AdapterRemoval v1.2-GG1 (Lindgreen, 2012). 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 (Gen-Bank accession JQ904290.1) using bwa (Li and Durbin, 2009). 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 (Picard, 2013 Felidae. Some phylogenetic studies of felids have been compromised by the inclusion of numt sequences in alignment (Davis et al., 2010). 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 Adap-terRemoval v1.2-GG1 (Lindgreen, 2012) 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.
We estimated the evolutionary relationships among pantherine cats using maximum likelihood in RAxML v8.0.14 (Stamatakis, 2014). 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 (Drummond et al., 2012). 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 (Xie et al., 2011). 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 (Drummond et al., 2006). We included relativerate 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 (Tseng et al., 2014), stem tiger Panthera zdanskyi (Mazak et al., 2011), and stem pantherine Panthera paleosinensis (Mazak, 2010). These fossils were used to specify uniform priors on the ages of corresponding nodes in the tree (Ho and Phillips, 2009). Previous estimates of the evolutionary timescale of cave lions were calibrated using fossil evidence from P. fossilis (Burger et al. 2004, Barnett et al. 2009), 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

Panthera uncia
Snow leopard NC010638 (Wei et al., 2011, Wei et al., 2009 Neofelis nebulosa Clouded leopard NC008450 (Wu et al., 2007) Felis sylvestris catus Domestic cat FCU20753 (Lopez et al., 1996) Felis sylvestris catus NUMT Nuclear pseudogene FCU20754 (Lopez et al., 1994, Lopez et al., 1996   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.

Radiocarbon dating
Sample YG401.410 We obtained two separate dating results for sample YG401.410. The dates were within three hundred 14 C 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 (Bocherens et al., 2011). 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 (2010) may be resolved with further sampling.

Sample F-2678/70
The hair sample has previously been radiocarbon dated by Kirillova et al. (2015). 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. (2015).

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 (Ma and Wang, 2014) and Asian lion (Bagatharia et al., 2013) produced a tree with P. spelaea as the sister taxon of P. l. persica, in contradiction with previous studies (Barnett et al., 2009, Burger et al., 2004, Ersmark et al., 2015. Analysis of the ATP8 gene, which has been characterised for both cytoplasmic and nuclear origin (Barnett et al., 2009) in Panthera cats, revealed the presence of a numt sequence within the mitochondrial genome published by Ma and Wang (2014). Given that this mitochondrial genome sequence was assembled from published nuclear genome data (Cho et al., 2013), it is likely to include a significant proportion of nuclearderived 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 (Johnson et al., 2006, Barnett et al., 2005. 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 (Burger et al., 2004, Barnett et al., 2009).

Evolutionary history of the lion-like cats
Despite their global range and continued dominance of ecosystems in Africa, and until recently in Asia, the lionlike 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 (Barry, 1987), with recognisably leonine fossils known from Olduvai II at 1.4-1.2 million years ago (Hemmer, 2011). 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) (Hemmer, 2011, Turner and Antón, 1997, Lewis et al., 2010, Sabol, 2011. 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 (Groiss, 1996, Herrington, 1986 and with the jaguar (Simpson, 1941, Christiansen andHarris, 2009), most have realised its close connection to the modern lion (Goldfuß, 1810, Dawkins et al., 1866, Vereshchagin, 1971, Turner, 1984, Kurtén, 1985, Sotnikova and Nikolskiy, 2006. Previous genetic studies (Burger et al., 2004 have used the first appearance of the ancestral cave lion (Panthera fossilis) Foronova, 2014, Peretto et al., 2015) 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) (Uphyrkina et al., 2001), or tiger (P. tigris) (Luo et al., 2004). 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 (Buckley-Beason et al., 2006).
Mitochondrial data and the considerable morphological differences (Sotnikova and Nikolskiy, 2006) 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.