A genomic exploration of the early evolution of extant cats and their sabre-toothed relatives

Background: The evolutionary relationships of Felidae during their Early–Middle Miocene radiation is contentious. Although the early common ancestors have been subsumed under the grade-group Pseudaelurus, this group is thought to be paraphyletic, including the early ancestors of both modern cats and extinct sabretooths. Methods: Here, we sequenced a draft nuclear genome of Smilodon populator, dated to 13,182 ± 90 cal BP, making this the oldest palaeogenome from South America to date, a region known to be problematic for ancient DNA preservation. We analysed this genome, together with genomes from other extinct and extant cats to investigate their phylogenetic relationships. Results: We confirm a deep divergence (~20.65 Ma) within sabre-toothed cats. Through the analysis of both simulated and empirical data, we show a lack of gene flow between Smilodon and contemporary Felidae. Conclusions: Given that some species traditionally assigned to Pseudaelurus originated in the Early Miocene ~20 Ma, this indicates that some species of Pseudaelurus may be younger than the lineages they purportedly gave rise to, further supporting the hypothesis that Pseudaelurus was paraphyletic.


Introduction
Based on the fossil record, the origin of Felidae, colloquially referred to as cats, is well resolved to the Middle-Late Oligocene (~30-27 million years ago (Ma)) 1 .However, the subsequent radiation of Felidae in the Early-Middle Miocene (~23-15 Ma) is less well understood, encompassing the evolution and radiation from the common ancestor of a number of species subsumed under the paraphyletic grade-group Pseudaelurus, which was widely distributed across Europe, Asia, and North America 1 .The putative paraphyly of Pseudaelurus is exemplified by the fact that it likely includes not only the early ancestors of modern cats (Felinae) but also those of the extinct sabretooths (Machairodontinae).
In contrast to the early radiation of Felidae, the phylogenetic relationships among extant cats are relatively well understood; extant cats are believed to share a most recent common ancestor in the Late Miocene (~11 Ma) 2 .However, despite this, there is uncertainty surrounding the constituents of the long-stem lineage of Felinae after its split from Machairodontinae (which extends through the Early and Middle Miocene) until its radiation in the Late Miocene.The evolutionary relationships within Machairodontinae are even less well understood, and increased knowledge of the evolutionary relationships both between Felinae and Machairodontinae, as well as within Machairodontinae itself, would provide important insights to help resolve the complex early radiation of Felidae.
The last surviving members of the Machairodontinae belonged to the genus Smilodon (tribe Smilodontini).While once widespread across the continents of North America (S.fatalis) and South America (S.populator), the genus went extinct ~10 thousand years ago (kya) 3,4 .Smilodon are not known further north than southernmost Canada (~42 °N), and their range extended until the tip of the South American continent (~53 °S).
Apart from a few anomalous mass occurrences in tar pits (e.g.Rancho La Brea, USA and Talara, Peru), Smilodon is relatively rare in the fossil record, which is not uncommon for an apex carnivore.Moreover, when it is present, there are only limited remains 5 .
The evolutionary relationships of Smilodon to other extinct and extant large cats are not resolved.An early ancient DNA study using both mitochondrial and nuclear DNA placed Smilodon within Felinae 6 .However, this was later shown to likely reflect contaminant DNA from a domestic cat 7 .A later study using complete mitochondrial genomes found Smilodon to be highly divergent from both Felinae (with an estimated divergence of ~20 Ma), and another extinct Machairodontinae lineage (Homotherium (tribe Homotheriini)), diverging ~18 Ma 8 .However, conclusions based exclusively on a single maternally inherited locus can be biased by interspecific hybridisation and incomplete lineage sorting, which has been well-documented in living cats 9 .
To date, only three studies present authentic DNA from Smilodon 8,10,11 , all of which is mitochondrial DNA.Although there is a large collection from Rancho La Brea, retrieving endogenous DNA from the material is not considered feasible 7 , greatly restricting the number of specimens available for DNA study.The relative absence of studies of Smilodon likely in part reflects their rarity in the fossil record outside of tar pits.In addition, as Smilodon are assumed to have been mixed wood-edge ambush predators 12 , the majority of fossils are found in open-air sites with poor DNA preservation; cave sites (which afford better preservation) are rare.Finally, the detrimental effects of temperature at equatorial and near-equatorial latitudes further exacerbate DNA preservation.
To elucidate the evolutionary relationship of Smilodon to extant and extinct Felidae, and to gain insights into the early radiation of Felidae, we sequenced a draft nuclear genome of a single Smilodon populator individual from the Ultima Esperanza region of Chile.The specimen was radiocarbon dated to 13,182 ± 90 calibrated years before present (cal BP), and represents the oldest palaeogenome from South America, a region in which ancient DNA preservation is expected to be limited 13 .

Results and discussion
Using an African lion assembly as the reference genome, we successfully mapped a draft nuclear genome of a single Smilodon populator individual to an average genome-wide coverage of ~0.7x, with an average read depth of ~2x.We achieved this using a combination of both target enrichment and shotgun sequencing.We performed target enrichment using baits created from either DNA extracted from a modern lion (whole-genome capture), or based on the exome annotations of the lion reference genome (exome capture).The discrepancy between genomewide coverage and average read-depth likely reflects the use of captured data, and lack of a closely related reference genome.If we were to have a conspecific reference genome, we would expect a more even genome-wide coverage, more comparable to the read depth.

Amendments from Version 1
We have made some minor updates to the manuscript based on the reviewers' suggestions.This includes fixing some typos and adding two new tables to the main text from the extended data.

Any further responses from the reviewers can be found at the end of the article
We sequenced approximately 800 million reads from four independently constructed libraries (Table 1).Investigations into mapped reads showed high levels of fragmentation (average read length <90bp for all libraries) and high rates of C-T transitions across the reads (Figure S1, Extended data 14 ), both typical of authentic ancient DNA.Comparisons among libraries showed that the capture experiments did yield significant levels of endogenous DNA, but were not highly different from the shotgun data but had lower levels of complexity.Moreover, although the exome capture was able to increase the relative coverage of the exome, it also included a lot of non-coding whole genomic data (Table 1).
To investigate the topological placement of Smilodon, we computed a neighbour-joining (NJ) tree using genome-wide transversion pairwise distances, rooted using Crocuta crocuta (spotted hyena).Similar to the mitochondrial genome, we found Smilodon and Homotherium to be sister taxa, and Machairodontinae to be sister taxon to all extant cats (Figure S2, Extended data 14 ).Support for this topology was high; we repeated this analysis independently for the 100 longest scaffolds (27.2Mb -5.3Mb), and found identical placement of Smilodon for each scaffold.
We further tested the closer affinity of Smilodon and Homotherium to each other, relative to any living Felinae, using D-statistics 15 to test for topology.For this, we took advantage of the topological input requirement for the D-statistics test ([[[H1, H2], H3], Outgroup]).We computed the D-statistic with Smilodon and Homotherium as sister taxa (i.e. in the H1 and H2 positions) (Table S1, Extended data 14 ), each species of Felinae as H3, and the spotted hyena as Outgroup, and compared this to the D-statistic and significance from 0 (Z-score) when placing Smilodon and Homotherium paraphyletically (i.e. in the H1 and H3 positions) (Table S2, Extended data 14 ).We clearly see a much higher D-statistic (D= 0.54-0.67,Z= 286.6-442.5)when placing the sabre-toothed cats parapyletically, compared to when they are placed as sister taxa (D= 0.12-0.20,Z= 39.2-67.7).A high D-statistic could be interpreted as gene flow between either H1 and H3, or between H2 and H3.However, another possible explanation could be the incorrect input topology.Although paraphyly of the sabre-toothed cats followed by gene flow could explain the observed pattern, a more likely explanation for the higher D-statistic when placing the sabre-tooths paraphyletically, would be the monophyly of the sabre-tooths as seen in our NJ tree.This result further supports Smilodon and Homotherium as more closely related to one another, than either are to any extant cat species.
To further contextualise this relationship, we built a dated phylogenetic tree (Figure 1).However, due to the low coverage of our Smilodon genome, traditional methods for phylogeny dating are likely unsuitable.Therefore, we devised a method to overcome this by computing pairwise genetic drift distances using F-statistics (in our case F2) (Table S3, Extended data 14 ) 16 .The pairwise F2 statistics were built into an unrooted NJ tree (Figure S3, Extended data 14 ).The relationships recovered in this tree were the same as those obtained using the transversional genetic distances.
To calibrate our phylogenetic tree, we estimated the average rate of genetic drift between all pairs of genera within Felinae based on the F2 results and previously calculated divergence dates 17 .We made a number of assumptions about our data, including (i) a strict molecular clock, (ii) constant population sizes through time, and (iii) that Felinae drift rates are similar to those in Machairodontinae.
Using the F2 statistics and average drift rate, we estimated a divergence date between Machairodontinae and Felinae of ~22.1 Ma (Figure 1), similar to the estimate of ~22.5 Ma calculated using a high-coverage exome dataset 17 .Although the correlation between these two results was not unexpected, as our tree was calibrated using the average within Felinae drift calculated using the divergence times of Barnett et al. 17 , their similar divergence estimates suggest extrapolating Felinae drift rates to other subfamilies of Felidae is valid.
We further tested for the robustness of this method by both recalculating the divergence times within Felinae, and by downsampling two Felinae species (Caracal caracal -caracal and Prionailurus bengalensis -leopard cat) to ~0.7x.We recovered similar divergence estimates to those previously reported, and to those produced without downsampling, providing confidence in the use of this methodology for low-coverage genomes (Tables S3 and S4, Extended data 14 ).Using this methodology, we estimated the divergence date of Smilodon/Homotherium to be ~20.65 Ma (95% CI 26.07-15.25 Ma) (Figure 1; Table S5, Extended data 14 ).Furthermore, although based on lowcoverage data and a number of data assumptions, the congruence of our results with previous mitogenome-based estimates from the same individuals 8 , despite the use of different data and calibration methods, further adds confidence to our method of phylogenetic assessment.
Species traditionally assigned to Pseudaelurus originated in the Early Miocene, ~20 Ma.Given the deep divergence between Machairodontinae and Felinae, as well as within Machairodontinae, this indicates that some species of Pseudaelurus are younger than the lineages they purportedly gave rise to.This is strong support for the hypothesis that Pseudaelurus as a grade-group is paraphyletic.Furthermore, current phylogenies using well-known characters (such as length of upper canines and presence of serrations or crenulations) resolve the Smilodontini-Homotheriini divergence to only ~11-10 Ma, while older machairodonts constitute a stem lineage of uncertain relationship 18 .Therefore, given the problematic nature of machairodont phylogenetics, the new molecular information analysed here provides an important reference point for identifying morphologically Smilodontini or Homotheriini characters in specimens from >11 Ma, to help resolve machairodont phylogeny back to the early Miocene.
We next investigated the evolutionary relationships between Smilodon and Felinae by searching for signatures of gene flow, using two independent methods (f3 16 and D3 19 ).The results were assessed using simulated data (Tables S6, S7, and S8, Extended data 14 ).As lineages leading to Smilodon and Homotherium diverged relatively soon after the Machairodontinae/ Felinae split, we were able to test for ancient gene flow (up to 20 Ma) between Smilodon or Homotheirum and stem Felinae.
Ancient gene flow may have prevented the diverging lineages in early Felidae from accumulating obvious morphological differences, thus preventing the confident phylogenetic placement of lineages during this time period.
Moreover, we tested for more recent gene flow events between either Smilodon or Homotherium and the entire Panthera big cat lineage (all Panthera species grouped together), due to their potential size overlap, as well as with single species that may have more recently met Smilodon due to spatial proximity (i.e. from South America: Panthera onca (jaguar), Puma concolor (puma), and Leopardus pardalis (ocelot)).
We did not find any indication of gene flow between any of the lineages tested, regardless of method (Tables S6 and S7, Extended data 14 ).The lack of more recent gene flow events is not surprising, due to the relatively ancient divergence of Machairodontinae and Felinae.However, the lack of ancient gene flow signatures do not necessarily mean that gene flow was not present during the early divergence of these lineages.Rather, it may result from the inadequacy of the methods in uncovering such events.
We assessed the power of the D3 statistic to detect gene flow, which occurred at different times (20 Ma-50 kya), using simulated data.When using a simple demographic model of constant population size, mutation rate, and recombination rate on error-free simulated data, we detect significant levels of gene flow back to ~16 Ma (Table S8, Extended data 14 ).However, empirical data do not always behave in such a simple manner, and a variety of factors may influence the results of the D3 statistic.These include, but are not restricted to: violations of the infinite sites model, different mutation rates across lineages, ancestral populations structure, and introgression with unsampled lineages 19 .However, our results suggest that D3 may be suitable for highly divergent lineages with ancient gene flow events and therefore violations of the infinite sites model may not be as problematic.
Thus, although we are unable to exclude the possibility of ancient gene flow during the early radiation of Felidae, we are somewhat more confident that there was no more recent gene flow between either Smilodon or Homotherium and Felinae within the last 16 Ma years.If very ancient gene flow had occurred prior to 16 Ma, it is likely that recombination and genetic drift would have either highly fragmented the introgressed regions, or completely removed them from contemporary Felinae individuals.This would render the regions too small to be detected with current methods, which use genome-wide summary statistics or regions of phylogenetic incongruence with known evolutionary relationships, to infer gene flow.
Our study exemplifies how even a draft palaeogenome from an extinct species can provide important information into their evolutionary history.Through the sequencing of a single Smilodon populator genome, we provide insights into Felidae's early radiation in the Early-Middle Miocene (~23 -15 Ma), which could not be uncovered using genetic data from extant species alone.

Sample information
Specimen ZMA20.042 is from the Naturalis museum, Leiden, the Netherlands, and was radiocarbon dated to 11,335±30 uncalibrated years before present (2-sigma range of 13,269-13,095 calibrated years before present) (Stafford: UCIAMS-142836).We calibrated the radiocarbon date using Calib v7.04 20 using the int13.14ccalibration curve.It has been identified as a left tibia of Smilodon populator and is part of the Kruimel collection, an assortment of megafaunal remains recovered from the Última Esperanza region of Chile, most likely from the site of Cueva del Milodón.This and other specimens of Smilodon from the Kruimel Collection were described and figured by McDonald and Werdelin (2018); specimen ZMA20.042 is presented in Figure 4.6D 5 .We additionally included genomic information from another extinct sabre-toothed cat, 12 extant cat species, and the spotted hyena (Table 2).
Ancient DNA extraction, library preparation, and sequencing Samples of cortical bone were taken from the long bone element (approx. 1 cm 3 ) using a Dremel drill, and reduced to powder in a Mikrodismembrator.Two independent DNA extractions were performed as described in Orlando et al. 21in a dedicated ancient DNA laboratory, with negative controls.We built each DNA extract and negative control into genomic libraries using the NEB E6070 kit following a modified version of the protocol used by Vilstrup et al. 22 .Briefly, the extract (30 µl) was end-repaired and cleaned using a MinElute column, the collected flow-through was adapter-ligated and cleaned using a QiaQuick column, and the adapter fill-in reaction was performed on the flowthrough.).

Leopardus pardalis Ocelot
For each library, we performed two independent indexing PCR amplifications (Veriti thermal cycler, Applied Biosystems) in a 50 µl volume reaction, using 25 µl of library, 25 µl PCR master mix, and 12 cycles of PCR reactions.The final concentrations in the PCR master mix were 1.25 U Accu-Prime™ Pfx DNA Polymerase (Invitrogen, Cat # 12344-024), 1x AccuPrime™ Pfx reaction mix (Invitrogen, Cat # 12344-024), 0.4 mg/ml BSA, 120 nM primer InPE, 120 nM of a multiplexing indexing primer containing a unique six-nucleotide index code (Illumina -sequences TGCAGG, CGATGA, GCGAGA, or CAGCAC).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, and a final extension step at 68°C for 7 min.Indexed libraries were checked for presence of DNA on a 2% agarose gel before purification using the QIAquick column system (Qiagen, Cat # 28104) and quantification on an Agilent 2100 BioAnalyzer.Quantified libraries were communally pooled in equimolar ratios and sequenced as 100 bp single-end reads on an Illumina HiSeq2000 platform at the Danish National High-Throughput Sequencing Centre and 100 bp paired-end reads on an Illumina Hiseq2000 at BGI Copenhagen.

Genome capture
We assessed the shotgun-sequenced libraries for endogenous content and selected the libraries with the highest levels of endogenous DNA for two sets of capture experiments.The first set used biotinylated RNA probes transcribed from fresh DNA extracted from a modern lion, for the purpose of enriching the entire nuclear genome (whole-genome capture).The second method used biotinylated RNA baits, assembled based on the exonic annotations of lion genomic data 23 (exome capture).Both types of baits were generated by Arbor Biosciences (Ann Arbor, MI, USA) and carried out using the myBaits target enrichment kit and the instructions described in manual V3.
After capture and cleanup, enriched libraries were re-amplified for further sequencing using either a Phusion polymerase (New England Biolabs, Cat # M0530S) or a KAPA HiFi HotStart polymerase (Roche, Cat # KK2801 07959052001) with primers IS5_reamp.P5 and IS6_reamp.P7 over 14 cycles 24 .We quantified the resultant enriched libraries on a TapeStation 2200 instrument and sequenced them on an Illumina Hiseq2000 at the Danish National High-throughput Sequencing Centre.

Data processing pipeline
Post-sequencing read processing of the Smilodon populator was performed using the PALEOMIX v1.2.5 pipeline 25 .
Adaptor sequence removal and trimming of low-quality bases (BaseQ < 5 or Ns) was done with AdapterRemoval v2.0.0 26 .This step also removed all reads shorter than 30 bp in length or with more than 50 bp of missing data.Trimmed reads were mapped against an African lion reference genome (NCBI Accession JAAVKH000000000 27 ) with BWA-MEM v0.7.5a 28 , utilising default parameters.PCR duplicates were identified and filtered based on the 5'-end mapping coordinate using Picard v2.18.0.GATK v3.8.0 29 was used to perform an indel realignment step to adjust for increased error rates at the end of short reads in the presence of indels.In the absence of a curated dataset of indels, this step relied on a set of indels identified in the specific sample being processed.Read damage patterns were assessed and base quality scores recalibrated around read terminal damage patterns using mapDamage v2.0.5 30,31 .
Data processing of the extant Felidae species, excluding Puma concolor and Leopardus pardalis, followed the same pipeline with the following minor adjustments; no minimum read length cut-off or missing data cut-off was applied during the adapter trimming step, and bases were not recalibrated using mapDamage.These steps were removed as the data from the extant species would not display the highly fragmented and damage patterns found in ancient DNA.The Puma concolor and Leopardus pardalis samples had Illumina adapter and short sequences trimmed using skewer 0.2.2 32 but followed the same protocol as the other extant species for the rest of the processing steps.
We further checked for support of the genome-wide topology by building a distance matrix for the 100 longest scaffolds independently, resulting in 100 independent distance matrices based on 171,281 -3,344,878 sites.The resultant distance matrices were converted into NJ trees using fastME v2.1.5 34.
We converted the haploid call file into a PLINK file format using the haplo2plink command from the ANGSD toolsuite.
Using the resultant PLINK file, we calculated F2 statistics 15 for each pairwise combination of our Smilodon populator individual, Homotherium latidens, and 12 Felinae species using the popstats.pyscript 36 .We used genetic drift calculated via F-statistics as opposed to absolute pairwise distance as it should be more suitable for ancient DNA data 37 .From these pairwise F2 comparisons, we built a distance matrix that was converted into a newick tree file using PHYLIP v3.696 neighbor 38 for visualisation.
From the distance matrix, we computed the Felinae average rate of genetic drift (F2) to be 0.000305 per million years.For this, we calculated the average F2 between pairwise comparisons of all genera within Felinae that had divergence time estimates available in Barnett et al. 17 .These included Acinonyx, Caracal, Felis, Lynx, Neofelis, Panthera, and Prionailurus, and resulted in 21 comparisons (Table S3, Extended data 14 ).We used the formula of (F2 × 0.5)/divergence time to estimate the rate of genetic drift per million years for all 21 comparisons.We calculated the mean of these 21 comparisons, giving us the Felinae mean rate of genetic drift per million years.This mean rate was used in conjunction with the previously calculated average F2 to calibrate the tree including the Smilodon populator, Homotherium latidens, puma, and ocelot.
To test for the robustness of our method to low-coverage data, we downsampled both Caracal and Prionailurus to comparable coverage to Smilodon populator (~0.7x) using SAMtools v1.6 39 , and recomputed the F2 statistics.We used the same average rate of genetic drift to estimate the divergence of these genera from their closest relatives (i.e.Prionailurus -Felis, and Caracal -Acinonyx/Felis/Lynx).For this analysis we assumed a known species tree and divergences within Felinae previously found in Barnett et al. 17 , a genomewide constant mutation/drift rate, and no variation in drift rates between lineages.Uncertainty around the divergence between Smilodon and Homotherium was calculated using 1.96 x the standard deviation on either side of the F2 statistic.
The standard deviation was calculated using the standard error x √ N (where N = number of blocks using for jackknifing).

Assessing gene flow
We implemented two independent analyses to test for the presence of signs of gene flow between Machairodontinae and Felinae.We computed F3 statistics 16 to assess whether there were any signals of gene flow between a predefined triplet [[A,B],C] using the same PLINK file computed above for the F2 statistics.We used popstats.py 36for five independent triplet combinations, which we present as ((A,B),C) in the subsequent text.We placed Smilodon and Homotherium in the A and B positions, while alternating C.
First, we investigated signs of very ancient gene flow between Machairodontinae and Felinae by specifying all extant Felinae as one population.Next, we looked for gene flow with the Panthera spp.big cats by specifying all Panthera as a single population.Finally, we computed F3 three times independently to test for signs of very recent gene flow between Machairodontinae and either of three extant cat species occupying South America, which may have come into contact with Smilodon based on geography (Panthera onca (jaguar), Puma concolor (puma), and Leopardus pardalis (ocelot)).
To complement this analysis, we also computed D3 statistics 18 on the same five triplet comparisons.For this, we also produced a haplocall file in ANGSD using the same filtering parameters as above, but using a random base call (-dohaplocall 1), as opposed to the consensus base call, while specifying the same additional parameters specified above.The resultant haploid output was then converted into a geno file and run through the popgenWindows.pypython script.We ran the popgenWindows.pyscript using default parameters, specifying a window size of 1MB, and the minimum number of sites per window as 1kb.From this output, we calculated D3 for each window independently by applying the equation Mean values, standard deviations, and significance from 0 were measured in R v3.6.0 using the pnorm function 40 .

Evaluating D3 for detecting ancient gene flow
To evaluate the adequacy of the D3 method for detecting very ancient gene flow events (up to 20 Ma), we simulated three diploid sequences representing Homotherium (A), Smilodon (B), and a Felinae species (C), using msprime 41 .As input, we estimated the average transversion mutation rate within Felidae using the average pairwise distance between Homotherium and Felinae and divided this by two before multiplying it by the Machairodontinae/Felinae divergence time calculated above (22.1 Ma).We computed pairwise comparisons between Homotherium and each Felinae species included in the study using ANGSD and took the average pairwise distance from these comparisons.We excluded the Smilodon from this calculation due to the lower quality of the genome.
To calculate the pairwise distances we used the consensus identity by state parameter (-doIBS 2) in ANGSD and applied the following filters; -minmapq 30, -minq 30, -minind 15, -uniqueonly 1, -remove_bads 1, -domajorminor 1, -rmtrans 1, -minminor 0, -makematrix 1, and only included scaffolds >1MB in length.This gave us an average pairwise distance of 0.01207, which we converted into an average transversion mutation rate of 2.730769e-10 per year.We converted this to a generational mutation rate using a generation time calculated for lions of 6.5 years 42 ; if the generation time of the investigated species differed from that of lion, the mutation rate would adjust accordingly, and we therefore expect minimal impact of this parameter on the final results, especially as we ran the simulations specifying years as opposed to number of generations.For the recombination rate, we used a previously published sex-averaged recombination rate of 1.9 cM/Mb (1.9e-8) 43 .
Using the above information, we ran five independent simulations, each consisting of 2,000 1 Mb windows with a constant effective population size of 40,000 individuals, a generation time of 6.5 years, the above-mentioned mutation and recombination rates, a 5% pulse of migration (m=0.05), and a different timing of the migration pulse for each of the five runs (20 Ma, 18 Ma, 17 Ma, 16 Ma, 15 Ma, 10 Ma, 5 Ma, 50 kya).The scripts for these simulation runs can be found on GitHub.We calculated the D3 statistic from this output using the tskit toolkit 44 using the d3.py script.Significance was calculated as it was done for the empirical data above.
This project contains the following extended data: In the present work, the authors show the results of the sequencing of a draft nuclear genome of the largest sabre-toothed felid, Smilodon populator, with an age of around 13,000 years ago, making this the oldest palaeogenome from South America to date, a region known to be problematic for ancient DNA preservation.The analysis of this genome, together with genomes from other extinct and extant cats, allows the authors to suggest a series of very interesting phylogenetic proposals.Among them, they confirm a deep divergence (around 20.65 Ma) within Machairodontinae, something that was evident from the fossil record, but never stated from palaeogenomics.
I haven't found major concerns about this work.The methodology seems impeccable, the figure is informative, the taxonomic citations are correct, and the discussion is well supported by the data.The style of the text is clear, and the manuscript is correctly written.I just detected a grammatical mistake when citing S. populator on page 3, which is cited as "popular".Finally, I find the subject treated in the manuscript suitable for its indexing on Open Research Europe.

MINOR COMMENTS
Although fossils of Early and Middle Miocene felids are known since the XIX century, our understanding of their phylogenetic relationships has not changed greatly over this time.This was in part due to the fact that those fossils of primitive felids, mostly fragmentary dental and postcranial material, show a very similar morphology to that of extant species, and only sabretoothed taxa exhibit evident differences in their elongated and flattened upper canines or their huge carnassial teeth.In the present paper, Westbury et al. have sequenced the genome of one of the most iconic sabre-toothed felids, the giant Smilodon populator, and then they compared it to those of extant species of Felidae, as well as to that of another extinct sabre-toothed species, the relatively slender and cursorial Homotherium latidens.The results were expectable but very interesting: they not only show the strong divergence between Machairodontinae (sabre-toothed felids) and Felinae (extant felids), but also between the two analysed species of sabre-toothed felids, S. populator and H. latidens.These two derived forms were the last sabre-toothed felids that were present in the faunas, and both were extinct as recently as around 28,000-11,000 years.Smilodon and Homotherium illustrate two very different ways to be a sabre-toothed felid: S. populator, included in a lineage known as Smilodontini, was a strong-built animal, with short and powerful limbs and very elongated upper canines, whereas H. latidens was a more slender predator, with relatively long limbs (the hind limbs much shorter than the forelimbs) and shorter canines.Thus, the divergence found by Westbury et   The study confirms that Smilodon and Homotherium are sister lineages, previously suggested by mtDNA data, and, as the authors point out, these results have implications for resolving the phylogeny of the machairodonts based on morphological characters.Gene flow between lineages during the Early Miocene radiation of the Felidae was not detected, within machairodonts or between machairodonts and stem Felinae, potentially owing to available data and analytical approaches that were possible.An absence of gene flow between machairodonts and extant cats in the more recent past is not too surprising.
An interesting study illuminating the history of a charismatic group.Unfortunately, perhaps owing to genome coverage, the authors did not examine genomic adaptations as Barnett et al. 2020 did for Homotherium latidens.

Tables & Figures
It seems unnecessary to avoid any Tables in the main text.For example, having the table of species analysed (Table S10) in the Methods section would benefit clarity, as these data suddenly find mention in the data processing section, as well as being featured in multiple subsequent sections.
The mapping results (Table S1) are summarized in the Results & Discussion section, but the table could feature in the main text as well.
Plain language summary "assess" rather than "access"?

Introduction
The first sentence reads rather clunky.

Results & Discussion
The start of this section would benefit from more details about (1) baits used for capture (based on data available for the lion; exome capture missing here?), and (2) mapping approach used (i.e. to which species: lion).
Can their approach to detect gene flow be applied to some Felinae (e.g.Panthera) in their dataset to confirm that it does indeed function for more recent events?

Methods
Exome capture used DNA baits?(typo?) Why was the lion genome used as reference?Why not a genome with higher contiguity?Did this not potentially result in unnecessary loss of data during filtering (e.g.due to scaffold lengths)?
Is the work clearly and accurately presented and does it cite the current literature?Yes Is the study design appropriate and does the work have academic merit?Yes

Are sufficient details of methods and analysis provided to allow replication by others? Yes
If applicable, is the statistical analysis and its interpretation appropriate?Yes Are all the source data underlying the results available to ensure full reproducibility?Yes

Are the conclusions drawn adequately supported by the results? Yes
Competing Interests: No competing interests were disclosed.
I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

Figure 1 .
Figure 1.Dated Felidae phylogenetic tree based on genome-wide pairwise F2 statistics, calibrated using an average genetic drift rate estimated from within Felinae (yellow shading).Blue bar shows 95% confidence interval for the divergence between the Smilodon (red) Homotherium lineages, calculated from the F2 standard error.Smilodon illustration by Binia De Cahsan and included with permission.

Reviewer Expertise:
Palaeobiology, Systematics, Functional anatomy I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.Daniel Förster 1 Evolutionary Genetics Department, Leibniz Institute for Zoo and Wildlife Research (IZW), Berlin, Germany 2 Evolutionary Genetics Department, Leibniz Institute for Zoo and Wildlife Research (IZW), Berlin, Germany Westbury et al. use nuclear data of a Smilodon populator individual to resolve the phylogenetic relationship of this species with another extinct sabre-toothed cat (Homotherium latidens), and of the Machairodontinae (here, Smilodon and Homotherium) with the extant cats (Felinae).They illuminate the deeper divergences within this Carnivoran family, i.e. the Felidae.

Table 2 . List of the additional species included in this study
. Original sources and accession codes for all raw reads used in this study can be found in TableS9(see Extended data In this test, 'extant cat' was replaced with each living Felidae species included in this study.We used ANGSD v0.921 used to infer post-divergence gene flow, but it can also be caused by more recent common ancestry brought about by an incorrect predefined topology.Taking the latter into account, we performed a number of D-statistics tests including the topologies: [[[Smilodon, Homotherium], extant cat], Crocuta crocuta], [[[Smilodon, extant cat], Homotherium], Crocuta crocuta], and [[[Homotherium, extant cat], Smilodon], Crocuta crocuta].

Is the work clearly and accurately presented and does it cite the current literature? Yes Is the study design appropriate and does the work have academic merit? Yes Are sufficient details of methods and analysis provided to allow replication by others? Yes If applicable, is the statistical analysis and its interpretation appropriate? Yes Are all the source data underlying the results available to ensure full reproducibility? Yes Are the conclusions drawn adequately supported by the results? Yes
3l. fits well with the idea that the sabretoothed felids developed at least two different morphotypes very early in their evolution, a robust, forest-dweller type illustrated by the Middle Miocene Pseudaelurus quadridentatus and the Late Miocene Promegantereon ogygia, and a slenderer model, probably occupying more open environments, whose first well-known representative is the Late Miocene Machairodus aphanistus.These two lineages evolved in parallel, refining each of the morphotypes, occupying different habitats thanks to different locomotor adaptations, thus allowing the presence in the faunas of two large predators through habitat segregation.Concerning the early divergence between sabre-toothed felids and felines, the present work by Westbury et al. points towards a date of around 22.1-22.5Ma.This is interesting, as we have little data on the felids of that moment, but all the available fossils suggest that the split between felines and macairodontines occurred at that time.There are several European localities with fossils slightly younger than that age (around 18 Ma) that can be identified as early felines, as their dentitions exhibit short canines of a rounded section, although they also show some very primitive features such as the presence of m2; and besides these forms, there are other fossils that show an incipient development of the macairodontine model illustrated by the laterally flattened upper canines of Pseudaelurus quadridentatus.Now seems clear that this latter genus should be considered as the first known sabre-toothed felid, as it has laterally flattened upper canines, a relatively robust first metacarpal, and an overall postcranial anatomy very similar to that of the younger Promegantereon ogygia (Beaumont, 1975 1 ; Salesa et al., 2010 2 ).Nevertheless, as the present work by Westbury et al. remarks, the genus Pseudaelurus has traditionally included felids with evident differences in their dental morphology, both from North America and Eurasia; the paraphyly of this genus was already suggested byWerdelin et al. (2010)3but this new work reinforces this idea.And finally, it is very interesting that, although Westbury et al. remark that they were unable to exclude the possibility of ancient gene flow during the early radiation of Felidae, they are somewhat more confident that there was no more recent gene flow between either Smilodon or Homotherium and Felinae within the last 16 Ma.Felinae and Machairodontinae developed a body plan devoted to preying upon different kinds of prey, and their dental and skeletal morphology suggests an early independent evolution of each model.This absence of gene flow also suggests the existence of several ethological and physiological differences between felines and sabretoothed felids whose trace unfortunately vanished when the latter disappeared.3. Werdelin L, Yamaguchi N, Johnson WE, O'Brien SJ: Phylogeny and evolution of cats (Felidae).Biology and Conservation of Wild Felids -Oxford University Press.2010.59-82 Competing Interests: No competing interests were disclosed.