Evolutionary history of Carnivora (Mammalia, Laurasiatheria) inferred from mitochondrial genomes

Abstract

The order Carnivora, which currently includes 296 species classified into 16 families, is distributed across all continents. The phylogeny and the timing of diversification of members of the order are still a matter of debate. Here, complete mitochondrial genomes were analysed to reconstruct the phylogenetic relationships and to estimate divergence times among species of Carnivora. We assembled 51 new mitogenomes from 13 families, and aligned them with available mitogenomes by selecting only those showing more than 1% of nucleotide divergence and excluding those suspected to be of low-quality or from misidentified taxa. Our final alignment included 220 taxa representing 2,442 mitogenomes. Our analyses led to a robust resolution of suprafamilial and intrafamilial relationships. We identified 21 fossil calibration points to estimate a molecular timescale for carnivorans. According to our divergence time estimates, crown carnivorans appeared during or just after the Early Eocene Climatic Optimum; all major groups of Caniformia (Cynoidea/Arctoidea; Ursidae; Musteloidea/Pinnipedia) diverged from each other during the Eocene, while all major groups of Feliformia (Nandiniidae; Feloidea; Viverroidea) diversified more recently during the Oligocene, with a basal divergence of Nandinia at the Eocene/Oligocene transition; intrafamilial divergences occurred during the Miocene, except for the Procyonidae, as Potos separated from other genera during the Oligocene.

About the Authors:

Alexandre Hassanin

Roles Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Resources, Writing – original draft

* E-mail: [email protected]

Affiliation: Institut de Systématique, Évolution, Biodiversité (ISYEB), Sorbonne Université, MNHN, CNRS, EPHE, UA, Paris, France

ORCID logo https://orcid.org/0000-0002-4905-8540

Géraldine Veron

Roles Resources, Writing – review & editing

Affiliation: Institut de Systématique, Évolution, Biodiversité (ISYEB), Sorbonne Université, MNHN, CNRS, EPHE, UA, Paris, France

Anne Ropiquet

Roles Resources, Writing – review & editing

Affiliation: Faculty of Science and Technology, Department of Natural Sciences, Middlesex University, London, United Kingdom

Bettine Jansen van Vuuren

Roles Resources, Writing – review & editing

Affiliation: Department of Zoology, Centre for Ecological Genomics and Wildlife Conservation, University of Johannesburg, Johannesburg, South Africa

Alexis Lécu

Roles Resources, Writing – review & editing

Affiliation: Parc zoologique de Paris, Muséum national d’Histoire naturelle, Paris, France

Steven M. Goodman

Roles Resources, Writing – review & editing

Affiliation: Field Museum of Natural History, Chicago, IL, United States of America

Jibran Haider

Roles Resources, Writing – review & editing

Affiliations Institut de Systématique, Évolution, Biodiversité (ISYEB), Sorbonne Université, MNHN, CNRS, EPHE, UA, Paris, France, Department of Wildlife Management, Pir Mehr Ali Shah, Arid Agriculture University Rawalpindi, Rawalpindi, Pakistan, Forest Parks & Wildlife Department Gilgit-Baltistan, Skardu, Pakistan

Trung Thanh Nguyen

Roles Resources, Writing – review & editing

Affiliation: Institut de Systématique, Évolution, Biodiversité (ISYEB), Sorbonne Université, MNHN, CNRS, EPHE, UA, Paris, France

Introduction

The order Carnivora is composed of 296 extant species [1] currently ranged into two suborders: the Caniformia which includes nine families, namely the Canidae (dog-like species), Ailuridae (red panda), Mephitidae (skunks and stink badgers), Mustelidae (weasels, badgers, martens, otters, etc.), Odobenidae (walrus), Otariidae (eared seals), Phocidae (earless seal), Procyonidae (raccoons, coatis, kinkajous, etc.) and Ursidae (bears); and the Feliformia which is represented by seven families, namely the Felidae (cat-like species), Eupleridae (Malagasy carnivorans), Herpestidae (mongooses), Hyaenidae (hyenas), Nandiniidae (African palm civet), Prionodontidae (Asiatic linsangs) and Viverridae (civets and genets). The oldest known fossils of Carnivora have been found in the late Paleocene; they belong to the extinct families Miacidae and Viverravidae, and have small body-size, comparable to extant weasels and martens [2–4]. The timing of the emergence of the crown carnivorans and their relationships to Paleocene and Eocene fossils are still unresolved. However, there is consensus for supporting two higher fossil taxa: the Carnivoraformes, which are composed of the crown group plus the stem family Miacidae, which is probably paraphyletic; and the Carnivoramorpha, which groups the Carnivoraformes and the Viverravidae [4,5].

The phylogeny of several carnivoran families has been extensively studied based on mitochondrial and nuclear data, including the Felidae, the Mustelidae and the Ursidae [6–8], while other families remain poorly studied, specifically the Eupleridae, Herpestidae, Mephitidae, Procyonidae and Viverridae. The phylogeny and timescale of diversification of the Carnivora have been studied by Eizirik et al. [9] using a molecular supermatrix of 7,765 base pairs (bp) containing 14 nuclear genes for 50 species, which represents less than 17% of the species diversity of the order Carnivora. The evolutionary history of carnivorans has also been inferred using a supertree approach by Nyakatura and Bininda-Emonds [10]. These two studies were based on different methods and data, and different fossils were used as calibration points. Although several nodes showed similar ages in the two studies, such as the most recent common ancestor (MRCA) of the Carnivora (59.2 Mya versus 64.9 Mya), Arctoidea (42.6 Mya versus 47.8 Mya) and Pinnipedia (24.5 Mya versus 22.4 Mya), some nodes were highly discordant, including the MRCA of the Caniformia (48.2 Mya versus 61.2 Mya), Feliformia (44.5 Mya versus 53.2 Mya), Feloidea (Felidae + Prionodontidae) (33.3 Mya versus 52.9 Mya) and Canidae (7.8 Mya versus 16.3 Mya). Most other studies to date have focused on either the Caniformia [8] or the Feliformia [11] and the timing of diversification in some families appears highly elusive or uncertain in the absence of a molecular timescale based on a high diversity of species from all carnivoran families.

With the development of next-generation sequencing (NGS) technologies, the number of mitochondrial genomes available in the international nucleotide databases has considerably increased during the last decade [12]. For the order Carnivora, there are currently more than 2,400 complete mitogenomes, and some of these were sequenced from Pleistocene fossils, such as polar bear [13], giant short-faced bears [14], cave lion [15] or saber-toothed cats [16]. This notably large and diversified dataset offers an excellent opportunity to better understand the evolutionary history of the order Carnivora, as problematic sequences and taxonomic issues can be more easily detected, and more importantly, as many fossils can be included as calibration points for estimating divergence times. The latter aspect is particularly relevant considering the paleontological record of Carnivora has been significantly improved over the last 10 years, with the discovery of several key fossils [17–21].

Here we analysed complete mitochondrial genomes to reconstruct phylogenetic relationships among carnivorans and to estimate divergence times. We sequenced 43 mitogenomes using various methods (Sanger sequencing of PCR products, NGS of long PCR products or Illumina shotgun sequencing). We also assembled eight mitogenomes from Sequence Read Archive (SRA) data. The 51 new mitogenomes belong to several families, namely Canidae (2), Eupleridae (6), Felidae (4), Herpestidae (12), Hyaenidae (1), Mustelidae (10), Otariidae (2), Phocidae (1), Prionodontidae (1), Procyonidae (3), Ursidae (1), Viverridae (7) for Carnivora, as well as Tapiridae (1) for the order Perissodactyla, which was included as an outgroup. These new mitogenomes were compared to all mitogenomes available for Carnivora in public repositories. At the intraspecific level, we selected only mitogenomes that were separated by more than 1% of nucleotide divergence and excluded those suspected to be of low-quality or from misidentified taxa. Our final alignment includes 220 taxa, which represent 2,442 mitogenomes. Our main objective was to build one of the largest time-trees of Carnivora estimated using a selection of fossil calibration points in order to provide further insights into the evolution of this broadly distributed and morphologically diverse order of mammals.

Material and methods

DNA extraction, amplification, sequencing and mitogenome assembly

Total DNA was extracted from cells, muscle or skin samples using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany). Details on the 43 samples extracted for this study are given in S1 Appendix. The mitochondrial genomes were sequenced using one of the three following approaches: Sanger sequencing of about 20 overlapping PCR products (length between 700 and 2000 bp); NGS of five overlapping long PCR products of around 4–5 kb; and Illumina shotgun sequencing.

In the first approach, PCR amplifications were carried out as previously described [22] using the primers listed in S2 Appendix. The amplicons were then sequenced in both directions by Eurofins MWG Operon (Ebersberg, Germany). Genomes were assembled with electropherograms of overlapping amplicons using Sequencher 5.1 (Gene Codes Corporation, Ann Arbor, MI, USA).

In the second approach, five overlapping PCR products of around 4–5 kb were amplified as previously described [23], and were sequenced at the “Service de Systématique Moléculaire” (UMS CNRS 2700, MNHN, Paris, France) using a Ion Torrent Personal Genome Machine (Thermo Fisher Scientific, Waltham, MA, USA).

All PCR products generated for this study were overlapping in the 5’ and 3’ regions with at least two other PCR products. As a consequence, the authenticity of each sequence was validated by double checking after verification of the perfect nucleotide identity with the two other sequences overlapping in the 5’ and 3’ regions (length between 50 and 500 bp).

The third approach was based on Illumina shotgun sequencing. DNA samples were quantified with a Qubit® 2.0 Fluorometer using the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA). Libraries were prepared using the TruSeq® Nano DNA Library Prep kit (Illumina, San Diego, CA, USA) after pooling 150 ng of total DNA of 10 species belonging to distant taxonomic groups (i.e. different phyla, classes, orders or families). Libraries were sequenced at the “Institut du Cerveau et de la Moelle épinière” (Paris, France) using a NextSeq® 500 system and the NextSeq 500 High Output Kit v2 (300 cycles) (Illumina).

The NGS reads generated with either Ion Torrent or Illumina sequencers were assembled by baiting and iterative mapping approach on Geneious® 10.2.2 (Biomatters Ltd., Auckland, New Zealand) using available mitochondrial references, including cytochrome b, cytochrome c oxidase subunit I, 12S and 16S rRNA genes. The 43 new mitochondrial genomes generated for this study were annotated using MITOS [24] and deposited in GenBank under accession numbers MW257198-MW257240.

Eight mitochondrial genomes were assembled from SRA downloaded from NCBI for the following species: Bassaricyon neblina (SRX1097850), Bassariscus sumichrasti (SRX1099089), Helogale parvula (SRR7637809), Lontra canadensis (SRR10409165), Mirounga angustirostris (SRR10331586), Mungos mungo (SRR7704821), Prionodon linsang (ERR2391707) and Zalophus wollebaeki (SRR4431565).

Mitochondrial alignments

The 51 mitochondrial genomes assembled in this study were compared with other genomes available in the NCBI nucleotide database (see details in S1 Appendix) for the different families. At the intraspecific level, we selected only mitochondrial haplotypes separated by more than 1% of divergence. Our final alignment includes 218 mitogenomes of Carnivora representing a large taxonomic diversity in the following 16 families (number of mitogenomes in parentheses): Ailuridae (2), Canidae (21), Eupleridae (7), Felidae (46), Herpestidae (15), Hyaenidae (4), Mephitidae (3), Mustelidae (39), Nandiniidae (1), Odobenidae (1), Otariidae (12), Phocidae (30), Prionodontidae (2), Procyonidae (6), Ursidae (21) and Viverridae (19). A tapir (Tapirus terrestris) and a pangolin (Phataginus tricuspis) were used to root the carnivoran tree, as these two taxa represent different Laurasiatherian orders, Perissodactyla and Pholidota, respectively, which were shown to be closely related to Carnivora in previous molecular studies [25,26].

The 220 mitochondrial genomes were aligned using AliView 1.22 [27]. Ambiguous regions for primary homology were excluded from the alignment. For that reason, the control region was completely removed, as well as some parts of rRNA and tRNA genes and a few nucleotides at the 5’- and 3’-regions of some protein-coding genes. To limit the impact of missing data, we also removed from the alignment all indels (insertions or deletions) detected in only one genome. The final alignment, named mtDNA, contains 220 taxa and 14,892 bp, representing 89% of the reference mitogenome for Canis lupus (NC_002008). Two other datasets were derived from mtDNA: (1) mtDNA-Tv (transversions only), in which the nucleotide G was replaced by A and the nucleotide T by C; and (2) PCG-mtDNA (10,809 nt), in which all regions other than protein-coding genes were removed (i.e. 12S and 16S rRNA genes and tRNA genes), as well as the ND6 gene (because it is located on the opposite strand of other protein-coding genes). The three datasets used in this study are available at https://osf.io/cfx8r/.

Analysis of base composition

The alignment of the protein-coding genes of 220 mitochondrial genomes (PCG-mtDNA dataset) was used to calculate the frequency of the four bases (A, C, G and T) at each of the three codon-positions (S3 Appendix). The 12 variables measured were then summarized by a principal component analysis (PCA) using the FactoMineR package [28] in R version 3.5.3 (from http://www.R-project.org/). The strand bias in nucleotide composition was studied at third codon-positions of the PCG-mtDNA dataset by calculating the relative frequencies of A and T nucleotides (AT3 skew = [A—T] / [A + T]) and the relative frequencies of C and G nucleotides (CG3 skew = [C—G] / [C + G]) [29–31].

Phylogenetic analyses

Two datasets (mtDNA and mtDNA-Tv) were analysed with probabilistic methods for tree reconstruction using the resources available from the CIPRES Science Gateway [32]. The Bayesian analyses were done with MrBayes 3.2.7 [33] using the two following models: GTR+I+G for mtDNA and JC69+I+G for mtDNA-Tv. The posterior probabilities (PP) were calculated using 10,000,000 Metropolis-coupled MCMC generations, tree sampling every 1000 generations and a burn-in of 25%.

To examine the phylogenetic signal within the mtDNA dataset, we also performed Bayesian analyses (with the same parameters) on 10 half-overlapping sub-datasets (i–x) of the about the same length (i.e., 2978 or 2980 bp), corresponding to the following positions: (i) 1–2978; (ii) 1489–4466; (iii) 2979–5956; (iv) 4467–7444; (v) 5957–8934; (vi) 7445–10422; (vii) 8935–11912; (viii) 10423–13400; (ix) 11913–14892; and (x) 13401–14892 + 1–1488. The use of half-overlapping sub-datasets (sliding window of ≈ 2978 bp) implies that all nucleotide sites of the total mtDNA alignment are represented twice in these Bayesian analyses. The lists of bipartitions obtained from Bayesian analyses of the 10 sub-datasets were transformed into a weighted binary matrix for supertree construction using SuperTRI v57 [34]. Each binary character corresponds to a node, which was weighted according to its frequency of occurrence in one of the 10 lists of bipartitions. In this manner, the SuperTRI method takes into account both principal and secondary phylogenetic signals. The SuperTRI bootstrap percentages (SBP) were obtained from PAUP 4* version 4b10 [35] after 1000 bootstrap replicates of the MRP (Matrix Representation with Parsimony) matrix of 3,398 binary characters (reconstructed under SuperTRI v57 [34]). Here, the SuperTRI analyses were conducted to test for phylogenetic signal along the mtDNA genome. If a robust node in the Bayesian tree (PP ≥ 0.95) is recovered with high SBP (≥ 95%) and repeated in most of the 10 Bayesian trees reconstructed from the half-overlapping sub-datasets of the mtDNA dataset, this signifies that the phylogenetic signal is present all along the mtDNA genome. If a node in the Bayesian tree is recovered with low SBP (< 95%) and repeated in less than five of the 10 Bayesian trees reconstructed from the half-overlapping sub-datasets of the mtDNA dataset, this indicates that the phylogenetic signal is weak or confined to a few fragments of the mtDNA genome. If there is a robust topological conflict between Bayesian and SuperTRI results, this suggests that at least one of the studied genomes was partially contaminated by a mitochondrial DNA sequence from another species or by a nuclear DNA sequence of mitochondrial origin (Numt). An example has been previously reported for the mitochondrial genomes of domestic goat [36].

Molecular dating

Divergence times were estimated on the CIPRES Science Gateway [32] using the mtDNA dataset and the Bayesian approach implemented in BEAST v.2.4.7 [37]. Twenty-one fossil calibration points were selected for molecular dating (Table 1). Most of these were interpreted from the fossil record using maximum (Max) and minimum (Min) ages. We applied two strategies for fossil calibration: (1) a uniform distribution between Max and Min on the calibrated node ages; or (2) a log-normal distribution on the calibrated node ages using Min as offset, M = Max—Min / 4 and S = 0.926 (to match the 97.5% quantile to Max) (see details in Table 1). The second strategy relies on the fact that minimum ages are generally more accurate and reliable than maximum ages because younger fossils are more abundant and more accurately dated than older fossils as a consequence of taphonomic processes and dating methods [38]. We applied a GTR + I + G model for the mtDNA alignment (with a proportion of invariants of 0.425) and a relaxed-clock model with uncorrelated lognormal distribution for substitution rates. Node ages were estimated using a calibrated Yule speciation prior and 2.108 generations, with tree sampling every 2,000 generations and a burn-in of 25% generations. MCMC mixing efficiency and convergence were assessed using the ESS values (>200) in Tracer v.1.7.1. The chronogram was reconstructed with TreeAnnotator, which is included in the BEAST package [37].

[Figure omitted. See PDF.]

Table 1. Maximum (Max) and minimum (Min) age calibrations (in millions of years ago, Mya) used for molecular dating analyses (with either uniform or log-normal prior distributions [U/L] on calibrated node ages interpreted from the fossil record).

https://doi.org/10.1371/journal.pone.0240770.t001

Results

Variation in base composition

The base composition (frequency of the nucleotides A, C, G and T) was analysed at the three codon-positions of the PCG-mtDNA dataset (S3 Appendix). The 12 variables measured for 220 taxa were summarized by a principal component analysis (PCA) based on the first two principal components (PC), which contribute 45.80% and 24.48% of the total variance, respectively (Fig 1). The variables factor map shows that the variance can be explained by similar differences in base composition at the three codon-positions. Some taxa have a mtDNA genome characterized by a divergent base composition.

[Figure omitted. See PDF.]

Fig 1. Variation in base composition of the mitogenomes of Carnivora.

The PCG-mtDNA dataset was used to calculate the frequency of the four bases (A, C, G and T) at each of the three codon positions, and the 12 variables measured were then summarized by a principal component analysis (PCA). The main graph represents the individual factor map based on 220 taxa. The families of Carnivora are highlighted by different colours. The small circular graph at the top left represents the variables factor map.

https://doi.org/10.1371/journal.pone.0240770.g001

Most members of the Mustelidae are found near the middle of the graph, except Mellivora because its mtDNA genome contains a higher percentage of guanine (15.65% versus “mean in other Mustelidae” [MoM] = 13.18%) and a lower percentage of adenine (28.35% versus MoM = 30.95%). This trend is observed at each of the three codon positions, and is more marked at third positions (S3 Appendix): G1 = 22.19% versus MoM = 20.87%, G2 = 12.19% versus MoM = 11.84%, G3 = 12.58% versus MoM = 6.82%; A1 = 30.55% versus MoM = 31.95%, A2 = 19.13% versus MoM = 19.41%, A3 = 35.38% versus MoM = 41.51%.

Within Viverridae, the mitogenome of the otter civet (Cynogale bennettii) is characterized by a higher percentage of thymine (32.90% versus “mean in other Viverridae” [MoV] = 30.73%) (Fig 1). This trend is observed at each of the three codon positions, and is more marked at third positions: T1 = 25.47% versus MoV = 23.79%, T2 = 42.71% versus MoV = 42.17%, T3 = 30.52% versus MoV = 26.21%.

Within Ursidae, the giant panda (Ailuropoda melanoleuca) shows a different base composition, with more adenines (30.63% versus “mean in other Ursidae” [MoU] = 29.55%) and thymines (‎31.02% versus MoU = 28.98%) than other species of the family. This trend is mainly explained by differences at third codon-positions, in which the giant panda shows more adenines (40.91% versus MoU = 38.46%) and thymines (‎28.52% versus MoU = 23.38%).

Most members of the Felidae, including Smilodon populator (Machairodontinae), are found near the middle of the graph, except Homotherium latidens (Machairodontinae) for which the mtDNA genome shows a very atypical base composition. At third codon-positions, Homotherium is characterized by a lower percentage of guanine (2.80% versus “mean in other Felidae” [MoF] = 6.25%) and a higher percentage of thymine (25.41% versus MoF = 20.79%). At second codon-positions, its mitogenome contains a lower percentage of thymine (40.80% versus “mean in other Felidae” [MoF] = 41.89%) and a higher percentage of guanine (12.58% versus MoF = 11.86%). At first codon-positions, the mtDNA of Homotherium shows a lower percentage of adenine (30.08% versus “mean in other Felidae” [MoF] = 32.07%), a higher percentage of guanine (22.65% versus MoF = 20.92%) and a higher percentage of guanine (23.07% versus MoF = 22.01%).

Within Canidae, the two extant tribes can be distinguished based on their nucleotide composition, as the mitogenomes of the Vulpini (fox-like canines) show more guanines (≥ 13.89% versus ≤ 13.28%) and less adenines (‎≤ 29.53% versus ≥ 29.83%) than those of the Canini (dog-like canines). This trend is mainly explained by a strong bias observed at third codon-positions, in which the Vulpini taxa have more guanines (G3 ≥ 8.39% versus ‎≤ 7.30%) and less adenines (A3 ‎≤ 37.18% versus ≥ 37.80%) than the Canini.

At third codon-positions, the Eupleridae and Hyaenidae have higher percentages of cytosine than other Feliformia (mean C3: Eupleridae = 36.78%; Hyaenidae = 35.55%; other Feliformia = 30.65%).

Pairwise mtDNA distances

Uncorrected pairwise distances were calculated on PAUP4* [35] using the mtDNA dataset (S4 Appendix). Intraspecific similarity distances are generally lower than 2%. However, there are several exceptions involving individuals from different geographic regions, which are or could be assigned to different subspecies or even species (see discussion): the extinct Japanese river otter (Lutra lutra, LC050126) versus extant representatives of Lutra lutra (2.1%); Ursus arctos isabellinus versus other brown bear subspecies (2.3–2.4%); Ursus thibetanus laniger versus Ursus thibetanus mupinensis (2.3%); the three samples of Melogale moschata (2.2–2.9%); the three samples of Viverricula indica (2.5–3.1%); Canis lupus familiaris versus Canis lupus chanco (2.6%); the two samples of Paradoxurus hermaphroditus (3.0%); the two samples of Leopardus pardalis (2.9%); the two mtDNA lineages identified in Prionailurus bengalensis (3.6%); the two samples of Herpestes javanicus (4.5%), although BLAST searches in NCBI suggest that the sequence NC_006835 belongs in fact to Herpestes auropunctatus (see discussion); the two samples of Mungos mungo (6.3%), but we suggest that the sample SRR7704821 may belong in fact to Mungos gambianus (see discussion).

Interspecific distances are generally higher than 2% of similarity. However, several species show more similar mitogenomes: Martes martes and Martes zibellina (1.9%); Phoca largha and Phoca vitulina (1.9%); Mustela putorius versus Mustela eversmannii (1.2%) and Mustela nigripes (1.6%); Felis catus and Felis silvestris (0.7%); Zalophus californianus and Zalophus wollebaeki (0.5%); Urocyon cinereoargenteus and Urocyon littoralis (0.4%). Within Arctocephalus forsteri, there are two mtDNA lineages: 17 mitogenomes are similar to the mitogenome of Arctocephalus australis (MG023139) (1.3%), whereas the 28 other mitogenomes are more divergent (2.1%). Within Ursus arctos, we found six mtDNA lineages showing between 1.1% and 2.4% of nucleotide divergence, but brown bears sampled on the Alaskan ABC islands (northern portion of the Alexander Archipelago) have mitogenomes which are highly similar to those of Ursus maritimus (0.4–0.5%).

Phylogeny of Carnivora

The Bayesian tree of Fig 2 was reconstructed from the mtDNA alignment. The results show that most nodes of the tree were highly supported (PP ≥ 0.95) and were similarly recovered as monophyletic entities using two other methods (nodes highlighted with a filled back circle in Fig 2), i.e. in the Bayesian tree reconstructed using the mtDNA-Tv alignment and JC69+I+G model (S5 Appendix) and in the Bootstrap 50% majority-rule consensus tree reconstructed from the MRP matrix of the SuperTRI analysis (S6 Appendix). Many of these nodes show maximal support in the SuperTRI analysis, as they were found monophyletic in all 10 Bayesian trees reconstructed from the 10 half-overlapping sub-datasets of the mtDNA alignment. For all these nodes, the phylogenetic signal is therefore robust across all parts of the mtDNA alignment.

[Figure omitted. See PDF.]

Fig 2. Phylogeny of Carnivora based on mitogenomes.

The Bayesian tree was reconstructed using the mtDNA dataset (220 taxa and 14,892 bp) and GTR+I+G model. The two outgroup species are not shown. Species names follow the classification of the IUCN [1]; the taxa written in red highlight the taxonomic issues discussed in the main text. The accession numbers of the 42 mitogenomes of Carnivora specially sequenced for this study are indicated in red. The eight mitogenomes here assembled from SRA data are shown in green. The blue circle associated to Mustela nudipes MH464792 indicates that the mitogenome was originally misassigned to Viverra tangalunga. Fossil species are followed by the symbol “†”. For each terminal taxon, the number of similar mitogenome(s) found in GenBank (pairwise distance < 1%) is indicated after the accession number. Dash branches indicate nodes supported by posterior probability (PP) < 0.95. Black circles indicate nodes that are also monophyletic in the two following trees: SuperTRI bootstrap 50% majority-rule consensus tree; and Bayesian tree obtained from the analysis of the mtDNA-Tv dataset and JC69+I+G model. Grey circles show nodes that are not found to be monophyletic with one of the two methods detailed above. White circles indicate nodes that are not monophyletic in both mtDNA-Tv and SuperTRI bootstrap consensus trees. No information was provided for the nodes highly supported by the SuperTRI analyses, i.e. which were found monophyletic in all the 10 Bayesian trees reconstructed from the 10 half-overlapping sub-datasets of the mtDNA dataset. For nodes less supported by the SuperTRI analyses, the number of Bayesian trees (< 10) showing the nodes is indicated.

https://doi.org/10.1371/journal.pone.0240770.g002

Among the nodes supported by all three methods of tree reconstruction and at least seven of the 10 Bayesian SuperTRI analyses are the following taxa: order Carnivora; suborders Caniformia and Feliformia; infraorders Arctoidea, Cynoidea, Feloidea and Viverroidea; all the 14 families represented by at least two members in our analyses (i.e. all families except Odobenidae and Nandiniidae); the subfamilies Euplerinae, Felinae, Galidiinae, Genettinae, Guloninae, Helictidinae, Hemigalinae, Herpestinae, Ictonychinae, Lutrinae, Machairodontinae, Melinae, Monachinae, Mungotinae, Mustelinae, Paradoxurinae, Phocinae, Tremarctinae, Ursinae and Viverrinae; the genera Catopuma, Felis, Genetta, Herpestes, Leopardus, Lutra, Lynx, Martes, Meles, Mirounga, Monachus, Mungos, Panthera, Paradoxurus, Phoca, Prionailurus, Prionodon, Puma, Urocyon, Vulpes and Zalophus; and the species Ailurus fulgens, Canis lupus, Civettictis civetta, Leopardus pardalis, Lutra lutra, Melogale moschata, Mustela sibirica, Panthera onca, Panthera pardus, Panthera tigris, Panthera uncia, Paradoxurus hermaphroditus, Prionailurus planiceps, Procyon lotor, Ursus maritimus, Ursus spelaeus, Ursus thibetanus and Viverricula indica. The results support the non-monophyly of four genera: (1) Arctocephalus, because Arctocephalus pusillus is divergent from Neophoca, Phocarctos and other species of Arctocephalus; (2) Canis, because Canis adustus is more distantly related to other Canis species than to Cuon alpinus; (3) Mustela, because Mustela frenata is the sister-species of Neovison vison, whereas all other species of Mustela are enclosed into a robust clade; and (4) Ursus, because Helarctos malayanus is closely related to Ursus americanus and Ursus thibetanus, whereas Ursus arctos, Ursus maritimus and Ursus spelaeus are grouped together. In addition, three species are not found monophyletic: (1) Arctocephalus forsteri is paraphyletic with respect to Arctocephalus australis; (2) Prionailurus bengalensis is paraphyletic with respect to Prionailurus viverrinus; and (3) Ursus arctos is paraphyletic with respect to Ursus maritimus.

BEAST chronogram inferred from the mtDNA alignment

Molecular estimates of divergence times are shown in the chronogram provided in Fig 3. The ages were inferred with BEAST using the mtDNA alignment and the 22 calibration points detailed in Table 1, including 21 fossil calibration points and a single molecular calibration point (used for the MCRA of Ferae). Two analyses were performed using either a uniform or log-normal prior distribution on the calibrated node ages (named “U approach” and “L approach”, respectively). The two BEAST chronograms show the same tree topology (Fig 3), which is similar to the one reconstructed under MrBayes (Fig 2). As expected, the chronogram inferred with the “L approach” show more recent estimates of divergence times (values highlighted in blue in Fig 3) than the chronogram inferred with the “U approach”.

[Figure omitted. See PDF.]

Fig 3. A molecular timescale for carnivoran evolution.

The estimates of divergence time were calculated under BEAST v.2.4.7 using the GTR+I+G model on the mtDNA dataset. The asterisks show the 21 fossil calibration points used for molecular estimation (see Table 1 for details). The chronogram, mean ages (values in black) and associated 95% confidence intervals (grey bars) were inferred using a uniform prior distribution for fossil calibration points (“U approach”). For comparison, the values in blue are mean ages estimated using a log normal prior distribution for fossil calibration points (“L approach”; see main text for details and discussion). Species names follow the classification of the IUCN [1]; the name of the taxa written in red have been changed following our results which suggest these taxonomic changes (see the discussion for details).

https://doi.org/10.1371/journal.pone.0240770.g003

Using the Geologic Time Scale v. 5.0 [48] for the correspondence between estimated divergence times and geologic epochs, the results suggest that the crown Carnivora divided into Caniformia and Feliformia during the early Eocene at around 52.7–46.7 Mya. Subsequently, the Canidae separated from other families in the early/middle Eocene at around 48.0–42.5 Mya, followed by the Ursidae and Pinnipedia in the middle/late Eocene at around 43.4–38.3 Mya and 41.4–36.5 Mya, respectively. The diversification of the Musteloidea began near the Eocene/Oligocene boundary at around 37.4–33.0 Mya, but the Mustelidae diverged from the Procyonidae in the early Oligocene at around 34.4–30 Mya. The family Phocidae separated from the Otarioidea in the Oligocene at around 27.2–22.9 Mya, whereas the separation between the Otariidae and Odobenidae occurred in the early Miocene at around 20.0–16.7 Mya. Within the Feliformia, the Nandiniidae diverged from other families at the Eocene-Oligocene transition at around 34.4–31.1 Mya. Then, the Feloidea and Viverroidea split in the Oligocene at around 29.7–27.0 Mya. The separation between the Felidae and Prionodontidae occurred in the Oligocene at around 27.4–24.7 Mya. Within the Viverroidea, the Viverridae diverged from other families in the Oligocene at around 28.2–25.8 Mya, followed by the Hyaenidae in the late Oligocene at around 24.5–22.9 Mya, whereas the split between the Eupleridae and Herpestidae took place in the Miocene at around 22.3–21.2 Mya. Within carnivoran families, the generic diversification occurred during the Miocene (all genera of Canidae, Herpestidae, Mephitidae, Procyonidae and Viverridae; most genera of Eupleridae, Felidae, Hyaenidae, Mustelidae, Phocidae and Ursidae) or more rarely during the Pliocene (most genera of Otariidae, excepting Callorhinus; Aonyx / Lutrogale; Caracal / Profelis; Histriophoca / Pagophilus; Hyaena / Parahyaena; Hydrurga / Leptonychotes; Galidictis / Mungotictis + Salanoia; Helarctos / Ursus) or Pleistocene (Halichoerus / Phoca / Pusa; Mungotictis / Salanoia).

Discussion

Mitogenomic variations at the species level

In 88% of the species for which at least two mitogenomes were available (65 out of 74), intraspecific distances were found to be less than 2%. In nine species, we detected mtDNA distances greater than 2% (S4 Appendix). Such a high mitogenomic divergence suggests that the samples may belong to two or more distinct species, because of either imperfect taxonomy or species misidentification (human error). Other explanations are however possible, such as mtDNA introgression from another species or high levels of mtDNA divergence due to strong female philopatry.

In mammals, females are generally philopatric, spending their lives close to their birthplace, whereas males typically undertake longer-distance dispersals during their lives [49–51]. Indeed, because they have to take care of their young, females tend to stay in areas where they can predict risks, as well as the resources. With time, this behaviour can result in a high spatial genetic structure of mtDNA variation since the mitogenome is transmitted maternally. The impact of female philopatry on the mtDNA evolution is expected to be more important if the preferred habitat-type is delimited by strong geographic barriers representing higher risks associated with dispersal of females (and their young), for example, open areas for small forest fruit bats [52] or large rivers and mountains for giraffes [53]. Female philopatry may therefore result in the selection of highly divergent mtDNA haplotypes in geographically isolated maternal lineages, whereas gene flow can be still maintained through long-distance male dispersals. Female philopatry can be advanced to explain the high intraspecific mtDNA distances found between populations of bears of the species Ursus arctos (six lineages, 1.1–2.4%) and Ursus thibetanus (four linages, 1.0–2.3%). This hypothesis is supported by previous phylogeographic studies showing strong geographic structure of mtDNA variation [54,55] and by field data on brown bears indicating that natal dispersal distances are five times smaller for females [56].

In several other species, we found two or three mitogenomic haplotypes or haplogroups which are separated by more than 2%. In the Eurasian otter (Lutra lutra), Waku et al. [57] sequenced three mtDNA genomes from extinct river otters of Japan: the two mtDNA genomes of otters from Honshu are similar to those found in extant Eurasian otters (< 1%), whereas the genome of the otter from Shikoku Island differs by 2.1%, suggesting that it belongs to a distinct subspecies (Lutra lutra nippon) or species (Lutra nippon). Since female philopatry versus male biased dispersal was well attested by both radiotracking and genetic data in European otters [58], we consider that nuclear data are needed to solve this taxonomic issue.

The three mtDNA haplotypes available for the small-toothed ferret badger (Melogale moschata) show comparatively high nucleotide distances, i.e., between 2.2% and 2.9%. The genus Melogale includes five species [1]: Melogale cucphuongensis known from a few animals collected in northern Vietnam [59,60], Melogale everetti on Borneo, Melogale moschata which is widely distributed from North-east India to Taiwan through China and Indochina, Melogale orientalis on Java, and Melogale personata, which is found in mainland Southeast Asia. In agreement with the recent study of Rozhnov et al. [61], our phylogenetic analysis based on cytochrome b (CYB) sequences (S7 Appendix) shows that Melogale cucphuongensis is the sister-group of a clade composed of three divergent mtDNA lineages: (i) Melogale personata, as represented by 16 samples from Vietnam; (ii) a first mtDNA haplogroup of Melogale moschata, which contains five samples from Vietnam (including our MW257240 sequence) and the reference mitogenome (NC_020644, unknown origin); and (iii) a second mtDNA haplogroup of Melogale moschata, which contains a sample from Taiwan (subspecies subaurantiaca) and a sample from China (MN400429, from Fei Zhou). Since the holotype of Melogale moschata was collected in the Guangdong Province of South China [62], our CYB tree suggests that the two samples from China and Taiwan represent the species Melogale moschata, whereas the samples from Vietnam currently assigned to Melogale moschata may belong to a new species, confirming that Vietnam is a key hot spot for the genus Melogale.

Three mitogenomic haplotypes were available for the small Indian civet (Viverricula indica), representing the three main mtDNA lineages previously identified [63], i.e. the subspecies V. i. indica from India, Madagascar (where it was introduced), Sri Lanka, West China and northern Indochina; V. i. rasse from southern Indochina to Java; and V. i. pallida, from East China to Taiwan. The mitogenomic distances are between 2.5% and 3.1%, suggesting that they may be treated as different species, but nuclear markers are needed to further investigate this taxonomic issue. Gaubert et al. [63] concluded that these three mitogenomic lineages diverged from each other in the Pliocene, between 3.2 and 2.6 Mya, but our dating estimates rather support that the divergences occurred in the Pleistocene, between 2.2–1.9 Mya and 1.6–1.4 Mya.

In the common palm civet (Paradoxurus hermaphroditus), the two mtDNA haplotypes differ by 2.9%, suggesting that they may come from two separate species. Phylogeographic studies based on both mitochondrial and nuclear sequences (CYB, control region, intron 7 of FGB) [64,65] supported the existence of three distinct species of Paradoxurus: Paradoxurus hermaphroditus (Indian and Indochinese regions), Paradoxurus musangus (mainland Southeast Asia, Sumatra, Java and other small Indonesian islands) and Paradoxurus philippinensis (Mentawai Islands, Borneo and the Philippines). Only the first two of these were represented in our mitogenomic study. Since they have overlapping distributions in mainland Southeast Asia, where it can be assumed that Paradoxurus hermaphroditus dispersed from the north of India while Paradoxurus musangus arrived from the south (see Fig 2 in [64]), the possibility of interspecific hybridization should be investigated in this region.

The two mitogenomes sequenced for the ocelot (Leopardus pardalis) differ by 3.0%, suggesting that this taxon may be split into two separate species, as proposed by Nascimento [66] based on a morphological analysis of 591 specimens of Leopardus spp. Using an alignment of the 5’-part of the mitochondrial control region, Eizirik et al. [67] have suggested a further division into four geographic groups. However, there is a 100-bp stretch of missing data in all their sequences, which is located in a region where we detected three motifs of 80 bp repeated in tandem in the mitogenome of Leopardus pardalis mitis and only two motifs in that of Leopardus pardalis pardalis. These repetitive sequences, named RS2 by Hoelzel et al. [68], were found in variable number in all families of Feliformia, from two in Prionodon pardicolor [69] to five in Nandinia binotata [70]. These variations in RS2 tandem repeats may pose serious problems of homology for DNA alignment. If RS2 repeats are removed from the alignment of Eizirik et al. [67], there is no robust signal for phylogenetic relationships within Leopardus pardalis (S8 Appendix). We recommend therefore to study the phylogeography of ocelot lineages with mitochondrial protein-coding genes, such as CYB and CO1 genes, or full mitogenomes. For taxonomy purpose, the analyses should be completed with nuclear markers, as Wultsch et al. [71] have shown evidence for female philopatry versus male-biased dispersal in ocelots sampled in Belize (Central America).

Two divergent mtDNA haplogroups were found for the leopard cat Prionailurus bengalensis (3.6%), corresponding to the Asian mainland leopard cat and Sunda leopard cat [72,73]. A species-level separation was confirmed by Y chromosome and whole-genome SNP studies [7,72]. Therefore, we follow Kitchener et al. [74] in recognizing two species of leopard cats: the mainland leopard cat Prionailurus bengalensis for Asian mainland leopard cats and the Sunda leopard cat Prionailurus javanensis for leopard cats from Java, Sumatra and Borneo. Previous molecular estimates have provided different ages for their divergence, between 2.67 and 0.92 Mya [7,72,73]. Although based on different calibration points, our estimate of 2.7–2.4 Mya corroborates that of Luo et al. [72]. This suggests that maternal linages were first isolated at the Pliocene/Pleistocene boundary, when the glacial/interglacial oscillations commenced. However, the mitogenomic tree does not support a sister-group relationship between the two species of leopard cats, as Prionailurus bengalensis is closer to Prionailurus viverrinus, the fishing cat, than to Prionailurus javanensis. This result contrasts with nuDNA trees showing a sister relationship of Prionailurus bengalensis and Prionailurus javanensis [7,72]. In agreement with the scenario proposed by Li et al. [7], such a discordance suggests that a mtDNA introgression occurred from Prionailurus bengalensis to Prionailurus viverrinus in the Pleistocene epoch, here dated at 2.1–1.8 Mya.

The two mitogenomes of the Javan mongoose (Urva javanica, previously named Herpestes javanicus, see [75] for details) differ by 4.5% in our alignment. BLAST searches in NCBI show that the reference genome sequenced for Urva javanica (NC_006835) belongs in fact to its sister-species, Urva auropunctata (small Indian mongoose), previously named Herpestes javanicus or Herpestes auropunctatus [75]. In agreement with this, the sample was collected in Fiji where the small Indian mongoose is known to have been introduced [76].

Two highly divergent mtDNA haplotypes (6.3%) were detected in banded mongooses (Mungos mungo) from Halle Zoological Garden (Germany) and from the San Diego Zoo Institute for Conservation Research (USA). The geographic origins of these animals are unknown, but such a high pairwise distance suggests that the two mitogenomes may represent the two species currently described in the genus Mungos: the banded mongoose (Mungos mungo), which is widely distributed in sub-Saharan Africa; and the Gambian mongoose (Mungos gambianus), which is endemic to West Africa. BLAST searches in NCBI show that our MW257205 sequence shares 99.9% of CYB identity with a mongoose from Kenya (AY928674) [77], confirming that it belongs to the species Mungos mungo. By contrast, the SRR7704821 sequence shows only 93% of identity with Mungos mungo AY928674, suggesting that it may rather belong to Mungos gambianus. Morphologically, the two species are similar but Mungos mungo is characterized by 10 to 15 dark dorsal stripes that are absent in Mungos gambianus. These marked differences render the hypothesis of misidentification unlikely. As an alternative, a mitochondrial introgression from Mungos gambianus to Mungos mungo may have occurred either in captivity or in the wild, as the two species can be found in sympatry in West Africa [1].

Shallow mitochondrial phylogeny

The mtDNA distances calculated between closely-related species are generally higher than 2% (73% of the comparisons detailed in S4 Appendix). However, there are several exceptions in the genera Arctocephalus, Felis, Martes, Mustela, Phoca, Urocyon, Ursus and Zalophus. Four main hypotheses can be advanced to explain the existence of similar mtDNA genomes in two putative species: species misidentification (human error), imperfect taxonomy (species synonymy), mtDNA introgression or recent speciation event during the Pleistocene [22].

Within Arctocephalus, two divergent mtDNA lineages (2.1%) were sequenced for the New Zealand fur seal (Arctocephalus forsteri) [78]. One of them was found to be more similar to the single mitogenome available for the South American fur seal, Arctocephalus australis (1.3%). The mitochondrial paraphyly of Arctocephalus forsteri suggests a mtDNA introgression from Arctocephalus australis at 0.8–0.7 Mya. Nuclear data should be sequenced to confirm this hypothesis and to further explore a possible taxonomic issue between South American and New Zealand fur seals.

Within Felis, the mitogenome of the domestic cat (Felis catus) is very similar to those of the wild cat, Felis silvestris (0.7%), including the African wild cat (Felis silvestris lybica), European wild cat (Felis silvestris silvestris) and Chinese desert cat (Felis silvestris bieti). These three subspecies and the domestic cat are either considered separate species [74] or subspecies of Felis silvestris [79,80]. Such a genome similarity between the domestic cat and wild cat was expected because the cat was domesticated from the Near Eastern wild cat Felis silvestris lybica [79] and can hybridize with wild forms [81]. Our results confirm that these taxa should be considered subspecies rather than species.

Within Martes, the pine marten (Martes martes) from western Europe and the sable (Martes zibellina) from Siberia and adjacent areas have quite similar mitogenomes (1.9% divergence), and our estimation of divergence time is 1.1–1.0 Mya, which is in agreement with Law et al. [8]. However, high levels of reciprocal mtDNA introgression have been detected in populations of northern Urals [82], where the two species are found in sympatry. In this hybrid zone, levels of gene flow should be further studied using nuclear markers to determine if introgression between the two taxa is asymmetric or not.

Within Mustela, several species show low levels of mitogenomic divergence (1.2-1-6%), including the western polecat (Mustela putorius) distributed in western Europe, the steppe polecat (Mustela eversmanii) found in eastern Europe and Asia, and the black-footed ferret (Mustela nigripes) endemic to North America. In addition, the mitogenome of the European mink, Mustela lutreola (MT304869, not included in our phylogenetic analyses) shares between 99.1% and 99.2% of identity with the four mitogenomes available for Mustela putorius. Our estimations indicate that these four closely related species of Mustela have diverged from each other after 1.1–1.0 Mya. However, Mustela putorius, Mustela eversmanii and Mustela lutreola can be found in sympatry in Europe [1], suggesting possible mtDNA introgression during the Pleistocene epoch. Cabria et al. [83] have shown evidence for mtDNA introgression from Mustela lutreola to Mustela putorius, resulting apparently from asymmetric interspecific hybridization, i.e. involving only females of European mink and males of polecat. Since the western and steppe polecats are known to occasionally hybridize in Eastern Europe [84], gene flow should be also studied between the two taxa.

Within Phoca, harbor seals (Phoca vitulina) and spotted seals (Phoca largha) are parapatric sibling species showing 1.9% of differences in our mitogenomic alignment. Although the two species are known to hybridize in captivity [1], they were found reciprocally monophyletic in previous mtDNA studies [85,86] and there is no evidence of mixed ancestry in wild populations found in sympatry (e.g. in Alaska [86]). The hypothesis of a recent Pleistocene speciation, at 1.1–1.0 Mya according to our estimation, is therefore the most likely explanation.

Within Urocyon, the island fox (Urocyon littoralis), restricted to six Channel Islands located off the coast of southern California (USA), and the mainland grey fox (Urocyon cinereoargenteus) share very similar mtDNA genomes (0.4%). In addition, the species Urocyon cinereoargenteus was found to be paraphyletic in the mitochondrial tree of Hofman et al. [87] due to the inclusive placement of Urocyon littoralis. Even if the island fox is approximately 25% smaller than the mainland grey fox [88], we suggest that the island fox should be rather treated as a subspecies of Urocyon cinereoargenteus. Nuclear genome comparisons are however needed to definitively clarify this taxonomic issue.

Within Ursus, bears living on the Alaskan ABC islands have mitogenomes which are very similar (0.4–0.5%) to those sequenced for extant polar bears and the fossil of Svalbard dated between 130 ka and 110 kya [13]. Two possible scenarios of mitochondrial introgression have been previously proposed: Hailer et al. [89] suggested that the ancestor of polar bears was introgressed by brown bears between 166 and 111 kya; whereas Hassanin [90] suggested that different populations of brown bears were introgressed by polar bears at two glacial periods of the Pleistocene, at 340 ± 10 ka in western Europe, and at 155 ± 5 ka on the ABC islands, and probably also in Beringia and Ireland based on ancient mtDNA sequences.

Within Zalophus, two species are currently recognized: the Californian sea lion (Zalophus californianus), which is found on the Pacific coasts of North America; and the Galápagos sea lion (Zalophus wollebaeki), which has been considered as a subspecies of Zalophus californianus by some authors [91]. Despite the low mitogenomic distance measured here between the two species (0.5%), Wolf et al. [92] have previously concluded that they are separate species because they were found reciprocally monophyletic with D-loop and CYB sequences and numerous private alleles were detected for both taxa at most of the 25 investigated microsatellite loci. Our molecular estimate of divergence time between these two species is 0.3–0.2 Mya, which is much more recent than previous estimations based on D-loop sequences, i.e. 0.8 Mya [93] and 2.3 ± 0.5 Mya [92].

Mitogenomic variations at the genus level

In mammals, most events of interspecific diversification at the genus level are more recent than 10 Mya and the great majority of them occurred during the Pliocene and Pleistocene epochs. This trend, previously reported in various genera of Cetartiodactyla [22] and Primates [94], is here confirmed in the following 25 carnivoran genera: Arctocephalus, Catopuma, Felis, Genetta, Leopardus, Lutra, Lynx, Martes, Meles, Melogale, Mirounga, Monachus, Mungos, Mustela, Panthera, Paradoxurus, Phoca, Prionailurus, Prionodon, Puma, Pusa, Ursus, Urva and Vulpes. Our analyses suggest however that two mustelid species separated from their congeneric representatives in the Middle Miocene: Martes pennanti at 14.2–12.4 Mya; and Mustela frenata at 13.4–11.8 Mya. In addition, the phylogenetic positions of these two species result in the polyphyly of the two genera Martes and Mustela, as previously found using both mtDNA and nuDNA markers [8]. Taken together, these results clearly indicate that Martes pennanti should be placed in its own genus Pekania, whereas Mustela frenata should be placed in the genus Grammogale. Indeed, the genus Grammogale was used to unite Mustela africana and Mustela felipei in a separate genus [95]. In addition, several molecular studies have shown that these two species fall into a robust clade with Mustela frenata and that they are the sister-group of Neovison vison [8,96]. We recommend therefore to include the four species Mustela africana, Mustela felipei, Mustela frenata and Neovison vison into the genus Grammogale.

Within the tribe Canini, our mitogenomic tree shows that Cuon alpinus (dhole) is the sister-group of all Canis species, excepting Canis adustus (side-striped jackal) which diverged from them at 8.7–7.4 Mya. The African wild dog (Lycaon pictus) occupies a more basal position, but this placement is only supported by 6/10 overlapping datasets in our SuperTRI analyses, indicating that an alternative position cannot be excluded. The mitochondrial phylogeny differs from published nuclear phylogenies [97,98] in which Canis adustus was found to be related to Canis mesomelas (black-backed jackal), but more divergent from the clade composed of other Canis species, Cuon and Lycaon. Following these results, in order to keep the genus Canis monophyletic, the species Canis adustus and Canis mesomelas should be placed in a different genus, which is Lupulella according to Viranta et al. [99]. The mito-nuclear discordance for the monophyly of the genus Lupulella was not discussed in previous studies. Our interpretation is that a mitochondrial introgression occurred at 6.2–5.2 Mya from an ancestor of Canis species to the lineage leading to Lupulella mesomelas. In agreement with the genomic study of Gopalakrishnan et al. [100], which concluded to pervasive gene flow among Canis species, our results suggest that interspecies hybridization has been also frequent in the early evolutionary history of canid genera in Africa.

Our chronogram in Fig 3 shows that no intergeneric divergence occurred during the Pleistocene epoch, except the split between the Malagasy euplerid genera Mungotictis and Salanoia, and the separation between the seal genera Halichoerus, Phoca and Pusa. A Pleistocene diversification was also found in previous molecular dating analyses of the Eupleridae and Phocidae [8,101]. In the absence of striking morphological feature for diagnosing these genera, we recommend synonymizing Halichoerus and Pusa with Phoca, as proposed in previous studies [102], and Mungotictis with Salanoia.

Changes in base composition and their impact on evolutionary rates

Mitogenomic rearrangements did not occur during the evolutionary history of the Carnivora. This deeply contrasts with other animal groups for which inversions of mitochondrial protein-coding genes and control region have resulted in many local or full reversals of base composition [30,31]. Despite this, five carnivoran species in our Bayesian tree (Fig 2) have a longer terminal branch than other representatives of the same family or superfamily: Cynogale bennetti within the Viverridae, Homotherium latidens within the Felidae, Mellivora capensis within the Mustelidae, Nasua nasua within the Procyonidae and Odobenus rosmarus within the Otarioidea. In addition, two higher taxa of Caniformia show longer branches: Otarioidea and Musteloidea. Since branch length is proportional to the number of substitutions, these long branches indicate either an acceleration in the rates of substitution or alternatively an important change in the pattern of substitution [103,104]. As shown in Fig 1, Cynogale, Homotherium and Mellivora show a divergent base composition. For Cynogale and Mellivora, we can therefore assume that a change in the pattern of substitution has taken place during their recent evolution. The case of Homotherium is more problematic. The mitogenome of this extinct genus was sequenced from a bone (YG 439.38) dated at >47,500 years [16]. Since ancient DNA molecules are known to exhibit a high rate of cytosine deamination, the sequences can contain artefactual C = >T and G = >A substitutions [105]. By comparing the mitogenomes of all felid species, this is exactly the pattern observed for third-codon positions, as Homotherium is characterized by less cytosines (28.94% versus 28.96–35.21% in other felid species), more thymines (25.41% versus 17.44–23.51%), less guanines (2.80% versus 4.69–8.41%) and more adenines (42.85% versus 38.41–42.79%). As a consequence, we conclude that sequencing errors introduced by DNA damage are the cause of the long branch of Homotherium in Fig 2.

The mitogenomes of Nasua and Odobenus have a base composition that is not so divergent from their closest relatives, suggesting that an acceleration of substitution rates is the cause of their long branch in Fig 2. As pointed out by Bromham [106], mammals evolving with faster evolutionary rates have generally a small body-size, but the cause of this body-size effect continues to be debated, because there are many possible mechanisms including shorter generation times, shorter lifespan, higher fecundity, larger litter size and higher metabolic rates. From this point, the walrus (Odobenus rosmarus) is intriguing because it is one of largest pinnipeds. Its field metabolic rate is 381.2 MJ/day, which is six to 32 times more than in other species of Pinnipedia [107]. Whereas most other pinnipeds hunt pelagic organisms, such as fish and cephalopod species, in the water column, the walrus feeds on benthic species and prefers molluscs, especially clams [108]. This special diet and associated behaviour, with long exposure in the cold coastal waters of the Arctic Ocean, may therefore explain the increased rate of mtDNA evolution detected in the walrus. Repetitive bottlenecks generated by climatic oscillations during the Pleistocene may be proposed as an alternative hypothesis. Nuclear genome comparisons between pinniped species will help to decipher between these two main hypotheses, as we can expect accelerated rates of substitution in both mtDNA and nuDNA genomes in case of repetitive bottlenecks.

Robust and reliable mitogenomic phylogeny for deep relationships within Carnivora

Previous molecular studies have shown that inferring deep phylogenetic relationships using mitogenomes can be problematic. An important issue concerns the reversals of strand mutational bias: detected in many animal phyla, such as Arthropoda, Chordata and Mollusca, it may be particularly misleading for phylogenetic reconstruction [30,31]. The mitogenome of mammals is not affected by this kind of bias, as its structure is highly conserved among the 27 mammalian orders. However, mutational saturation is a different issue: since the mtDNA evolves with higher mutation rates than the nuclear genome [109], it is more prone to multiple substitutions at the same site which leads, with time, to the disappearance of the phylogenetic signal [110,111]. Although mutational saturation may be particularly problematic for reconstructing Precambrian, Paleozoic and Mesozoic divergences, it is expected to have less impact for inferring more recent divergences, such as those during the Cenozoic. This explains why mitogenomic sequences have been largely used for inferring interfamilial and intrafamilial relationships in most mammalian groups, including Cetartiodactyla [22], Feliformia [11], Primates [94,112], Pholidota [113], Xenarthra [114], etc. In agreement with this view, all families and inter-familial levels shown in our mtDNA tree of Fig 2 received high support with all the three phylogenetic approaches (Bayesian trees of mtDNA and mtDNA-Tv datasets and SuperTRI analysis) indicating considerable stability of the topology. Interfamilial relationships are conformed to the nuclear tree published by Eizirik et al. [9], except the position of the Mephitidae as the sister-group of the Ailuridae. However, this relationship was also found by Law et al. [8] based on 46 genes (4 mitochondrial and 42 nuclear) and 75 species of Musteloidea.

Diversification of the Carnivora during the Cenozoic

The mtDNA phylogeny reconstructed here presents a good opportunity to study the evolution of the order Carnivora. Indeed, our dense taxonomic sampling allow us to include many fossils as calibration points for estimating divergence times. This point is particularly relevant as the fossil record of Carnivora has been significantly improved over the last 10 years, with the discovery of several key fossils [17–21]. As pointed out in Warnock et al. [115], the most effective means of establishing the quality of fossil-based calibrations is through a priori evaluation of the intrinsic palaeontological, stratigraphic, geochronological and phylogenetic data. We identified therefore 21 fossil calibration points (Table 1), including one tip calibration (the mitogenome of a late Pleistocene fossil of Ursus maritimus [13]) and 20 well-constrained calibration points having both minimum and maximum age boundaries interpreted from fossils with known position with respect to extant taxa included in our study. Molecular dating analyses were performed using either a uniform (U) or log-normal (L) prior distribution on the calibrated node ages. The “L approach” used here considers that minimum ages are generally more accurate and reliable than maximum ages because younger fossils are generally more abundant and precisely dated than older fossils as a consequence of taphonomy and dating methods [38]. As expected, the divergence times estimated with the “L approach” are more recent than those estimated with the “U approach” (Figs 3 and 4). The values calculated with the two approaches are given for all the nodes described below.

[Figure omitted. See PDF.]

Fig 4. Comparison with published chronograms on Carnivora.

The mean divergence times were here estimated with two approaches for the prior distribution on the calibrated node ages: (1) a uniform distribution between maximum and minimum boundaries (“U approach”, blue histograms); and (2) a log-normal distribution (“L approach”, green histograms) (see main text for details). The results were compared with mean ages inferred in Eizirik et al. [9] (red histograms) and Nyakatura and Bininda-Emonds [10] (orange histograms).

https://doi.org/10.1371/journal.pone.0240770.g004

The age estimated for the MRCA of Carnivora was 52.7 Mya with the “U approach” and 46.7 Mya with the “L approach”, and both are younger than previous estimates published by Eizirik et al. [9] (59.2 Mya) and Nyakatura and Bininda-Emonds [10] (64.9 Mya). However, our dating estimates fit well with the end of the warmest period of the Cenozoic era, from the Palaeocene–Eocene Thermal Maximum (PETM) at 56 Mya to the Early Eocene Climatic Optimum (EECO) at around 53–50 Mya [116–118]. In the fossil record, this period was marked by the first appearances of Primates, Perissodactyla and Cetartiodactyla in North America and Europe, as well as Carnivoraforms, a group formed by the crown group Carnivora plus the stem family ‘Miacidae’ [2,4,119]. During the early Eocene greenhouse world, rainforests spread on Earth from pole to pole [120,121], suggesting that early carnivorans were small arboreal species. Biogeographically, Carnivoraforms diversified in the three Laurasian continents during the early Eocene, with a possible origin in the late Paleocene of Asia [4]. Crown carnivorans originated in one of the three continents of the Northern Hemisphere. To know which one, new fossils need to be discovered in the middle Eocene to fill the gap between the Carnivoraforms of the early Eocene and the oldest carnivorans of the late Eocene [4].

The divergence times estimated for basal relationships in the Caniformia are similar to those published by Eizirik et al. [9] but much younger than those published by Nyakatura and Bininda-Emonds [10]: Caniformia = 48.0–42.5 Mya versus 48.2 Mya and 61.2 Mya respectively; Arctoidea = 43.4–38.3 Mya versus 42.6 Mya and 47.8 Mya respectively. For more recent nodes, our estimates are similar to the two previously cited studies: Musteloidea = 37.4–33.0 Mya versus 32.0 Mya and 35.7 Mya respectively; Mustelidae+Procyonidae = 34.4–30.3 Mya versus 27.4 Mya and 31.6 Mya respectively; and Pinnipedia = 27.2–22.9 Mya versus 24.5 Mya and 22.4 Mya respectively. By contrast, the ages inferred for the families are generally much younger in Eizirik et al. [9] (Fig 4), a result potentially explained by a poor species sampling (two species of Ursidae, three species of Canidae, etc.) and the low phylogenetic signal of nuclear protein-coding genes for the most recent nodes.

The ages estimated for basal divergences in the Feliformia are much younger than in Eizirik et al. [9] and Nyakatura and Bininda-Emonds [10]: Feliformia = 34.4–31.1 Mya versus 44.5 Mya and 53.2 Mya respectively; Feloidea = 27.4–24.7 Mya versus 33.3 Mya and 52.9 Mya respectively; and Viverroidea = 28.2–25.8 Mya versus 37.4 Mya and 31.8 Mya respectively. These results may be explained by the use of different calibration points: Feloidea = 34–20 Mya versus >28.5 Mya in Eizirik et al. [9]; Felidae = 20–14 Mya versus >31.15 Mya in Nyakatura and Bininda-Emonds [10]; and Viverridae = 34–14 Mya versus >25.72 Mya in Nyakatura and Bininda-Emonds [10]. The minimum age constraints used by Nyakatura and Bininda-Emonds [10] for Felidae and Viverridae seem however problematic because they were not extracted from the fossil record, but from the supertree of Bininda-Emonds et al. [122]. By contrast, the minimum boundary of 28.5 Mya used for Feloidea in Eizirik et al. [9] was interpreted from early Oligocene fossils assumed to be the oldest felids according to McKenna and Bell [123]. In Meredith et al. [26], who used a similar minimum boundary for Feloidea (> 28.3 Mya), the divergence times estimated for Feliformia were found to be highly variable between their eight molecular dating analyses, with mean ages between 35.81 and 43.07 Mya. Three early Oligocene fossils were used as calibration points in this study: Proailurus and Stenogale, which were assumed to be stem Felidae; and Palaeoprionodon, which was supposed to be a stem Prionodontidae. However, the phylogenetic positions of these fossils remain uncertain. Indeed, the three fossil genera formed the sister-group of Feloidea + Viverroidea + Herpestides antiquus† in the phylogenetic analyses of Solé et al. [3,4]. In the classification of Hunt [124], Palaeoprionodon was included in the subfamily Prionodontinae with Prionodon and other extant viverrid genera, such as Genetta and Pioana. Although Prionodon was treated as a member of the Viverridae by Hunt [124] and previous authors, it was then found to be the sister-group of the Felidae by Gaubert and Veron [125] which was confirmed by all more recent studies [9,10,26] (see also Fig 2). In Nyakatura and Bininda-Emonds [10], Paleoprionodon was used as a calibration point for the Viverridae (based on Hunt and Tedford [126]), although it is suggested to be either close to Prionodon or a stem Feliformia [3,124,127].

Our dating estimates suggests that the basal split between Nandinia and other genera of Feliformia occurred at the Eocene/Oligocene transition (34 Mya), when a brutal and global cooling of 5°C resulted in the extinction of many taxa and the appearance of several modern mammal lineages [128]. Accordingly, Nandinia may be the descendant of forest-adapted ancestors, whereas the other lineage of Feliformia may have evolved in response to this important climatic change by adapting to more open vegetation.

Conclusion

Based on a large taxonomic sampling (220 taxa represented by 2,442 mitogenomes) and 21 fossil calibration points (most of them having both minimum and maximum age boundaries and a position based on recent fossil revisions and phylogenies), our study proposes a new time-scaled phylogeny for Carnivora and provides further insights into the evolutionary history of this order. The age estimates for the Carnivora and the two suborders Caniformia and Feliformia fit well with global changes corresponding to the appearance of other mammal lineages. Moreover, our phylogenetic results, although largely similar to other recent phylogenies, suggest several taxonomic changes that would need to be confirmed using nuclear data.

Supporting information

S1 Appendix. Origin of the sequences used in this study.

https://doi.org/10.1371/journal.pone.0240770.s001

(PDF)

S2 Appendix. PCR primers used in this study.

https://doi.org/10.1371/journal.pone.0240770.s002

(PDF)

S3 Appendix. Base composition.

https://doi.org/10.1371/journal.pone.0240770.s003

(PDF)

S4 Appendix. Distances.

https://doi.org/10.1371/journal.pone.0240770.s004

(PDF)

S5 Appendix. Bayesian tree reconstructed using the mtDNA-Tv dataset (220 taxa and 14,892 bp) and JC69+I+G model.

https://doi.org/10.1371/journal.pone.0240770.s005

(PDF)

S6 Appendix. Bootstrap 50% majority-rule consensus tree reconstructed from the MRP matrix of the SuperTRI analysis.

https://doi.org/10.1371/journal.pone.0240770.s006

(PDF)

S7 Appendix. Cytochrome b phylogeny of the genus Melogale.

https://doi.org/10.1371/journal.pone.0240770.s007

(PDF)

S8 Appendix. D-loop phylogeny of the species Leopardus pardalis.

https://doi.org/10.1371/journal.pone.0240770.s008

(PDF)

Acknowledgments

We thank our colleagues and their institutions for access to samples and for their help during this study: N. Chai and the staff of the Ménagerie du Jardin des Plantes, the staff of the Teuk Chhou zoo, H.E. Oknha Nhim Vanda, F. Catzeflis, R. Cornette, C. Denys, G. Dobigny, A. Delapré and the Mammal collection of the Muséum national d’Histoire naturelle, A. Délicat, J. Fuchs, P. Gaubert, C. Hatten, Y. Varelides, the Phongsaly Provincial Agriculture and Forestry Office and the EU Phongsaly Forest Conservation and Rural Development Project, S. Heard, A. Kitchener from the National Museums Scotland, J.-P. Hugot, S. Laidebeure and the staff of the Parc zoologique de Paris, E. Leroy, Van Cung Nguyen, P. Nicolau-Guillaumet, F. Njiokou, I. Parker M.-L. Patou, le Parc Zoologique de Montpellier, T. Petit and the staff of the zoo de La Palmyre, X. Pourrut, P. Taylor and the Durban Natural Science Museum, the staff of Southport Zoo, Do Tuoc, B. Patterson, W. Stanley and J. Phelps from the Field Museum of Natural History, V. Volobouev and L. Woolaver.

Citation: Hassanin A, Veron G, Ropiquet A, Jansen van Vuuren B, Lécu A, Goodman SM, et al. (2021) Evolutionary history of Carnivora (Mammalia, Laurasiatheria) inferred from mitochondrial genomes. PLoS ONE 16(2): e0240770. https://doi.org/10.1371/journal.pone.0240770

References

1. The IUCN Red List of Threatened Species. Version 2020–1. https://www.iucnredlist.org. Downloaded on 19 March 2020.

2. Spaulding M., Flynn JJ. Phylogeny of the Carnivoramorpha: the impact of postcranial characters. J Vertebr Paleontol. 2012; 10: 653–677.

3. Solé F, Smith R, Coillot T, de Bast E, Smith T. Dental and tarsal anatomy of ‘Miacis’ latouri and a phylogenetic analysis of the earliest carnivoraforms (Mammalia, Carnivoramorpha). J Vertebr Paleontol. 2014; 34: 1–21.

4. Solé F, Smith T, de Bast E, Codrea V, Gheerbrant E. New Carnivoraforms from the latest Paleocene of Europe and their bearing on the origin and radiation of Carnivoraformes (Carnivoramorpha, Mammalia). J Vertebr Paleontol. 2016; 36: e1082480

5. Flynn JJ, Finarelli JA, Spaulding M. Phylogeny of the Carnivora and Carnivoramorpha, and the use of the fossil record to enhance understanding of evolutionary transformations; In Goswami A, Friscia A, editors. Carnivoran evolution: new views on phylogeny, form, and function. Cambridge UK: Cambridge University Press; 2010. pp. 25–63.

6. Kutschera VE, Bidon T, Hailer F, Rodi JL, Fain SR, Janke A. Bears in a forest of gene trees: phylogenetic inference is complicated by incomplete lineage sorting and gene flow. Mol Biol Evol. 2014; 31: 2004–17. pmid:24903145

7. Li G, Davis BW, Eizirik E, Murphy WJ. Phylogenomic evidence for ancient hybridization in the genomes of living cats (Felidae). Genome Res. 2016; 26: 1–11. pmid:26518481

8. Law CJ, Slater GJ, Mehta RS. Lineage diversity and size disparity in Musteloidea: testing patterns of adaptive radiation using molecular and fossil-based methods. Syst Biol. 2018; 67: 127–144. pmid:28472434

9. Eizirik E, Murphy WJ, Koepfli KP, Johnson WE, Dragoo JW, Wayne RK, et al. Pattern and timing of diversification of the mammalian order Carnivora inferred from multiple nuclear gene sequences. Mol Phylogenet Evol. 2010; 56: 49–63. pmid:20138220

10. Nyakatura K, Bininda-Emonds OR. Updating the evolutionary history of Carnivora (Mammalia): a new species-level supertree complete with divergence time estimates. BMC Biol. 2012; 10: 12. pmid:22369503

11. Zhou Y, Wang SR, Ma JZ. Comprehensive species set revealing the phylogeny and biogeography of Feliformia (Mammalia, Carnivora) based on mitochondrial DNA. PLoS One. 2017; 12: e0174902. pmid:28358848

12. Coissac E, Hollingsworth PM, Lavergne S, Taberlet P. From barcodes to genomes: extending the concept of DNA barcoding. Mol Ecol. 2016; 25: 1423‐1428. pmid:26821259

13. Lindqvist C, Schuster SC, Sun Y, Talbot SL, Qi J, Ratan A, et al. Complete mitochondrial genome of a Pleistocene jawbone unveils the origin of polar bear. Proc Natl Acad Sci USA. 2010; 107: 5053–7. pmid:20194737

14. Mitchell KJ, Bray SC, Bover P, Soibelzon L, Schubert BW, Prevosti F, et al. Ancient mitochondrial DNA reveals convergent evolution of giant short-faced bears (Tremarctinae) in North and South America. Biol Lett. 2016; 12: 20160062. pmid:27095265

15. Barnett R, Zepeda Mendoza ML, Rodrigues Soares AE, Ho SYW, Zazula G, Yamaguchi N, et al. Mitogenomics of the extinct cave lion, Panthera spelaea (Goldfuss, 1810), resolve its position within the Panthera cats. Open Quat. 2016; 2: 1–11.

16. Paijmans J, Barnett R, Gilbert M, Zepeda-Mendoza ML, Reumer J, de Vos J, et al. Evolutionary history of saber-toothed cats based on ancient mitogenomics. Curr Biol. 2017; 27: 3330–3336.e5. pmid:29056454

17. Abella J, Alba DM, Robles JM, Valenciano A, Rotgers C, Carmona R, et al. Kretzoiarctos gen. nov., the oldest member of the giant panda clade. PLoS One. 2012; 7: e48985. pmid:23155439

18. Grohé C, Bonis L de, Chaimanee Y, Blondel C, Jaeger JJ. The oldest Asian Sivaonyx (Lutrinae, Mustelidae): a contribution to the evolutionary history of bunodont otters. Palaeontol Electron. 2013; 16: 29A.

19. Tseng ZJ, Wang X, Slater GJ, Takeuchi GT, Li Q, Liu J, et al. Himalayan fossils of the oldest known pantherine establish ancient origin of big cats. Proc R Soc B. 2014; 281: 20132686. pmid:24225466

20. Berta A, Churchill M, Boessenecker RW. The origin and evolutionary biology of Pinnipeds: seals, sea lions, and walruses. Annu Rev Earth Planet Sci. 2018; 46: 203–228.

21. Valenciano A, Abella J, Alba DM, Robles JM, Álvarez-Sierra MA, Morales J. New Early Miocene Material of Iberictis, the oldest member of the wolverine lineage (Carnivora, Mustelidae, Guloninae). J Mammal Evol. 2020; 27: 73–93.

22. Hassanin A, Delsuc F, Ropiquet A, Hammer C, Jansen van Vuuren B, Matthee C, et al. Pattern and timing of diversification of Cetartiodactyla (Mammalia, Laurasiatheria), as revealed by a comprehensive analysis of mitochondrial genomes. C R Biol. 2012; 335: 32–50. pmid:22226162

23. Hassanin A, Bonillo C, Tshikung D, Pongombo Shongo C, Pourrut X, Kadjo B, et al. Phylogeny of African fruit bats (Chiroptera, Pteropodidae) based on complete mitochondrial genomes. J Zool Syst Evol Res. 2020; 58: 1395–1410.

24. Bernt M, Donath A, Jühling F, Externbrink F, Florentz C, Fritzsch G, et al. MITOS: improved de novo metazoan mitochondrial genome annotation. Mol Phylogenet Evol. 2013; 69: 313–9. pmid:22982435

25. Murphy WJ, Eizirik E, O’Brien SJ, Madsen O, Scally M, Douady CJ, et al. Resolution of the early placental mammal radiation using Bayesian phylogenetics. Science. 2001; 294: 2348–2351. pmid:11743200

26. Meredith RW, Janečka JE, Gatesy J, Ryder OA, Fisher CA, Teeling EC, et al. Impacts of the Cretaceous terrestrial revolution and KPg extinction on mammal diversification. Science. 2011; 334: 521–524. pmid:21940861

27. Larsson A. AliView: a fast and lightweight alignment viewer and editor for large data sets. Bioinformatics. 2014; 30: 3276–3278. pmid:25095880

28. Lê S, Josse J, Husson F. FactoMineR: An R Package for multivariate analysis. J Stat Softw. 2008; 25: 1–18.

29. Lobry JR. Properties of a general model of DNA evolution under no-strand-bias conditions. J Mol Evol. 1995; 40: 326–330. pmid:7723059

30. Hassanin A, Léger N, Deutsch J. Evidence for multiple reversals of asymmetric mutational constraints during the evolution of the mitochondrial genome of Metazoa, and consequences for phylogenetic inferences. Syst Biol. 2005; 54: 277–298. pmid:16021696

31. Arabi J, Judson ML, Deharveng L, Lourenço WR, Cruaud C, Hassanin A. Nucleotide composition of CO1 sequences in Chelicerata (Arthropoda): detecting new mitogenomic rearrangements. J Mol Evol. 2012; 74: 81–95. pmid:22362465

32. Miller MA, Pfeiffer W, Schwartz T. Creating the CIPRES Science Gateway for inference of large phylogenetic trees. In: New Orleans, LA: Proceedings of the Gateway Computing Environments Workshop (GCE); 2010: pp. 1–8.

33. Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Höhna S, et al. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol. 2012; 61: 539–542. pmid:22357727

34. Ropiquet A, Li B, Hassanin A. SuperTRI: a new approach based on branch support analyses of multiple independent data sets for assessing reliability of phylogenetic inferences. C R Biol. 2009; 332: 832–847. pmid:19748458

35. Swofford DL. PAUP*. Phylogenetic Analysis Using Parsimony (*and Other Methods). Version 4. Massachusetts: Sinauer Associates, Sunderland; 2003.

36. Hassanin A, Bonillo C, Nguyen BX, Cruaud C. Comparisons between mitochondrial genomes of domestic goat (Capra hircus) reveal the presence of numts and multiple sequencing errors. Mitochondrial DNA. 2010; 21: 68–76. pmid:20540682

37. Bouckaert R, Heled J, Kühnert D, Vaughan T, Wu CH, Xie D, et al. BEAST2: a software platform for Bayesian evolutionary analysis. PLoS Comput Biol. 2014; 10: e1003537. pmid:24722319

38. Crees JJ, Collen B, Turvey ST. Bias, incompleteness and the ’known unknowns’ in the Holocene faunal record. Phil Trans R Soc B. 2019; 374: 20190216. pmid:31679489

39. Tedford RH, Wang X, Taylor BE. Phylogenetic systematics of the North American fossil Caninae (Carnivora: Canidae). Bull Am Mus Nat Hist. 2009; 325: 1–218.

40. Banyue W, Zhanxiang Q. Notes on Early Oligocene ursids (Carnivora, Mammalia) from Saint Jacques, Nei Mongol, China. Bull Am Mus Nat Hist. 2003; 279: 116–124.

41. Bonis L de, Gardin A, Blondel C. Carnivora from the early Oligocene of the ‘Phosphorites du Quercy’ in southwestern France. Geodiversitas. 2019; 41: 601–621.

42. Soibelzon LH, Tonni EP, Bond M. The fossil record of South American short-faced bears (Ursidae, Tremarctinae). J South Am Earth Sci. 2005; 20: 105–113.

43. Rybczynski N, Dawson MR, Tedford RH. A semi-aquatic Arctic mammalian carnivore from the Miocene epoch and origin of Pinnipedia. Nature. 2009; 458: 1021–4. pmid:19396145

44. Dewaele L, Lambert O, Louwye S. A critical revision of the fossil record, stratigraphy and diversity of the Neogene seal genus Monotherium (Carnivora, Phocidae). Royal Soc Open Sci. 2018; 5: 171669. pmid:29892365

45. Sato JJ, Wolsan M, Minami S, Hosoda T, Sinaga MH, Hiyama K, et al. Deciphering and dating the red panda’s ancestry and early adaptive radiation of Musteloidea. Mol Phylogenet Evol. 2009; 53: 907–22. pmid:19699810

46. Werdelin L, Yamaguchi N, Johnson WE, O’Brien SJ. Phylogeny and evolution of cats (Felidae). In: Macdonald DM, Loveridge A, editors. The biology and conservation of wild felids. Oxford: Oxford University Press; 2010: pp. 59–82.

47. Werdelin L, Peigné S. Carnivora. In: Werdelin L, Sanders WJ, editors. Cenozoic mammals of Africa. University of California Press; 2010: pp. 603–657.

48. Walker JD, Geissman JW, Bowring SA, Babcock LE, compilers. Geologic Time Scale v. 5.0: Geological Society of America, https://doi.org/10.1130/2018.CTS005R3C. 2018.

49. Greenwood PJ. Mating systems, philopatry and dispersal in birds and mammals. Animal Behav. 1980; 28: 1140–1162.

50. Dobson FS. The enduring question of sex-biased dispersal: Paul J. Greenwood’s (1980) seminal contribution. Animal Behaviour. 2013; 85: ‏ 299–304.

51. Li XY, Kokko H. Sex-biased dispersal: a review of the theory. Biol Rev Camb Phil Soc. 2019; 94: 721–736. pmid:30353655

52. Hassanin A, Khouider S, Gembu GC, Goodman SM, Kadjo B, Nesi N, et al. The comparative phylogeography of fruit bats of the tribe Scotonycterini (Chiroptera, Pteropodidae) reveals cryptic species diversity related to African Pleistocene forest refugia. C R Biol. 2015; 338: 197–211. pmid:25636226

53. Petzold A, Hassanin A. A comparative approach for species delimitation based on multiple methods of multi-locus DNA sequence analysis: A case study of the genus Giraffa (Mammalia, Cetartiodactyla). PLoS One. 2020; 15: e0217956. pmid:32053589

54. Hirata D, Mano T, Abramov AV, Baryshnikov GF, Kosintsev PA, Vorobiev AA, et al. Molecular phylogeography of the brown bear (Ursus arctos) in Northeastern Asia based on analyses of complete mitochondrial DNA sequences. Mol Biol Evol. 2013; 30: 1644–1652. pmid:23619144

55. Wu J, Kohno N, Mano S, Fukumoto Y, Tanabe H, Hasegawa M, et al. Phylogeographic and demographic analysis of the Asian Black Bear (Ursus thibetanus) based on mitochondrial DNA. PLoS One. 2015; 10: e0136398. pmid:26406587

56. Steyaert SMJG, Endrestøl A, Hackländer K, Swenson JE, Zedrosser A. The mating system of the brown bear Ursus arctos. Mammal Review. 2012; 42: 12–34.

57. Waku D, Segawa T, Yonezawa T, Akiyoshi A, Ishige T, Ueda M, et al. Evaluating the phylogenetic status of the extinct Japanese otter on the basis of mitochondrial genome analysis. PLoS One. 2016; 11: e0149341. pmid:26938434

58. Quaglietta L, Fonseca V.C., Hájková P, Mira A, Boitani L. Fine-scale population genetic structure and short-range sex-biased dispersal in a solitary carnivore, Lutra lutra. J Mammal. 2013; 94: 561–571.

59. Nadler T, Streicher U, Stefen C, Schwierz E, Roos C. A new species of ferret-badger, Genus Melogale, from Vietnam. Zoologische Garten. 2011; 80: 271–286.

60. Li S, Yu GH, Liu S, Jin CS. First record of the ferret-badger Melogale cucphuongensis Nadler et al., 2011 (Carnivora: Mustelidae), with description of a new subspecies in southeastern China. Zool Res. 2018; 40: 575–579.

61. Rozhnov VV, Korablev MP, Abramov AV. Systematics and distribution of Ferret Badgers Melogale (Mammalia, Mustelidae) in Vietnam: first genetic data. Dokl Biol Sci. 2019; 485: 47–51. pmid:31197594

62. Storz JF, Wozencraft WC Melogale moschata. Mammalian Species. 1999; 631: 1–4,

63. Gaubert P, Patel RP, Veron G, Goodman SM, Willsch M, Vasconcelos R, et al. Phylogeography of the Small Indian Civet and origin of introductions to Western Indian Ocean Islands. J Hered. 2017; 108: 270–279. pmid:27940474

64. Patou ML, Wilting A, Gaubert P, Esselstyn JA, Cruaud C, Jennings AP, et al. Evolutionary history of the Paradoxurus palm civets–a new model for Asian biogeography. J Biogeogr. 2010; 37: 2077–2097.

65. Veron G, Patou ML, Tóth M, Goonatilake M, Jennings AP. How many species of Paradoxurus civets are there? New insights from India and Sri Lanka. J Zoolog Syst Evol Res. 2015; 53: 161–174.

66. Nascimento FO do. Revisão taxonômica gênero do Leopardus Gray, 1842, PhD Thesis, Sao Paulo: University of Sao Paulo. 2010. https://doi.org/10.11606/T.41.2010.tde-09122010-104050

67. Eizirik E, Bonatto SL, Johnson WE, Crawshaw PG Jr, Vié JC, Brousset DM, et al. Phylogeographic patterns and evolution of the mitochondrial DNA control region in two Neotropical cats (Mammalia, Felidae). J Mol Evol. 1998; 47: 613–624. pmid:9797412

68. Hoelzel AR, Lopez JV, Dover GA, O’Brien SJ. Rapid evolution of a heteroplasmic repetitive sequence in the mitochondrial DNA control region of carnivores. J Mol Evol. 1994; 39: 191–199. pmid:7932782

69. Hassanin A. The complete mitochondrial genome of the African palm civet, Nandinia binotata, the only representative of the family Nandiniidae (Mammalia, Carnivora). Mitochondrial DNA A DNA Mapp Seq Anal. 2016; 27: 904–905. pmid:24937569

70. Hassanin A, Veron G. The complete mitochondrial genome of the Spotted Linsang, Prionodon pardicolor, the first representative from the family Prionodontidae (Mammalia, Carnivora). Mitochondrial DNA A DNA Mapp Seq Anal. 2016; 27: 912–913. pmid:24937573

71. Wultsch C, Waits LP, Kelly MJ. A Comparative analysis of genetic diversity and structure in Jaguars (Panthera onca), Pumas (Puma concolor), and Ocelots (Leopardus pardalis) in fragmented landscapes of a critical Mesoamerican linkage zone. PLoS One. 2016; 11: e0151043. pmid:26974968

72. Luo SJ, Zhang Y, Johnson WE, Miao L, Martelli P, Antunes A, et al. Sympatric Asian felid phylogeography reveals a major Indochinese-Sundaic divergence. Mol Ecol. 2014; 23: 2072–2092. pmid:24629132

73. Patel RP, Wutke S, Lenz D, Mukherjee S, Ramakrishnan U, Veron G, et al. Genetic structure and phylogeography of the Leopard Cat (Prionailurus bengalensis) inferred from mitochondrial genomes. J Hered 2017; 108: 349–360. pmid:28498987

74. Kitchener AC, Breitenmoser-Würsten C, Eizirik E, Gentry A, Werdelin L, Wilting A, et al. A revised taxonomy of the Felidae: the final report of the Cat Classification Task Force of the IUCN Cat Specialist Group. Cat News. 2017; S11: 47–48.

75. Patou ML, McLenachan PA, Morley CG, Couloux A, Jennings AP, Veron G. Molecular phylogeny of the Herpestidae (Mammalia, Carnivora) with a special emphasis on the Asian Herpestes. Mol Phylogenet Evol. 2009; 53: 69–80. pmid:19520178

76. Veron G, Jennings AP. Javan mongoose or small Indian mongoose–who is where? Mammal Biol. 2017; 87: 62–70.

77. Koepfli KP, Jenks SM, Eizirik E, Zahirpour T, Van Valkenburgh B, Wayne RK. Molecular systematics of the Hyaenidae: relationships of a relictual lineage resolved by a molecular supermatrix. Mol Phylogenet Evol. 2006; 38: 603–620. pmid:16503281

78. Emami-Khoyi A, Paterson AM, Hartley DA, Boren LJ, Cruickshank RH, Ross JG et al. Mitogenomics data reveal effective population size, historical bottlenecks, and the effects of hunting on New Zealand fur seals (Arctocephalus forsteri). Mitochondrial DNA A DNA Mapp Seq Anal. 2018; 29: 567–580. pmid:28539070

79. Driscoll CA, Menotti-Raymond M, Roca AL, Hupe K, Johnson WE, Geffen E, et al. The Near Eastern origin of cat domestication. Science. 2007; 317: 519–523. pmid:17600185

80. Mattucci F, Oliveira R, Lyons LA, Alves PC, Randi E. European wildcat populations are subdivided into five main biogeographic groups: consequences of Pleistocene climate changes or recent anthropogenic fragmentation? Ecol Evol. 2016; 6: 3–22. pmid:26811770

81. Randi E. Detecting hybridization between wild species and their domesticated relatives. Mol Ecol. 2008; 17: 285–293. pmid:18173502

82. Rozhnov VV, Meschersky IG, Pishchulina SL, Simakin LV. Genetic analysis of sable (Martes zibellina) and pine marten (M. martes) populations in sympatric part of distribution area in the northern Urals. Russ J Genet. 2010; 46: 488–492.

83. Cabria MT, Michaux JR, Gómez-Moliner BJ, Skumatov D, Maran T, Fournier P., et al. Bayesian analysis of hybridization and introgression between the endangered European mink (Mustela lutreola) and the polecat (Mustela putorius). Mol Ecol. 2011; 20: 1176–1190. do: pmid:21244536

84. Abramov AV, Puzachenko AY, Tumanov IL. Morphological differentiation of the skull in two closely-related mustelid species (Carnivora: Mustelidae). Zoological Studies. 2016; 55: 1. pmid:31966146

85. Nakagawa E, Kobayashi M, Suzuki M, Tsubota T. Genetic variation in the harbor seal (Phoca vitulina) and spotted seal (Phoca largha) around Hokkaido, Japan, based on mitochondrial cytochrome b sequences. Zool Sci. 2010; 27: 263–268. pmid:20192695

86. Cordes LS, O’Corry-Crowe G, Small RJ. Surreptitious sympatry: exploring the ecological and genetic separation of two sibling species. Eco Evol. 2017; 7: 1725–1736. pmid:28331583

87. Hofman CA, Rick TC, Hawkins MT, Funk WC, Ralls K, Boser CL, et al. Mitochondrial genomes suggest rapid evolution of dwarf California Channel Islands foxes (Urocyon littoralis). PLoS One. 2015; 10: e0118240. pmid:25714775

88. Schutz H, Polly PD, Krieger JD, Guralnick RP. Differential sexual dimorphism: size and shape in the cranium and pelvis of grey foxes (Urocyon). Biol. J. Linn. Soc. Lond. 2009; 96: 339–353.

89. Hailer F, Kutschera VE, Hallström BM, Klassert D, Fain SR, Leonard JA, et al. Nuclear genomic sequences reveal that polar bears are an old and distinct bear lineage. Science. 2012; 336: 344–347. pmid:22517859

90. Hassanin A. The role of Pleistocene glaciations in shaping the evolution of polar and brown bears. Evidence from a critical review of mitochondrial and nuclear genome analyses. C R Biol. 2015; 338: 494–501. pmid:26026577

91. Aurioles-Gamboa D, Hernández-Camacho J. Zalophus californianus. The IUCN Red List of Threatened Species 2015: e.T41666A45230310.

92. Wolf JB, Tautz D, Trillmich F. Galápagos and Californian sea lions are separate species: Genetic analysis of the genus Zalophus and its implications for conservation management. Front Zool. 2007; 4: 20. pmid:17868473

93. Schramm Y, Mesnick SL, de la Rosa J, Palacios DM, Lowry MS, Aurioles-Gamboa D, et al. Phylogeography of California and Galápagos sea lions and population structure within the California sea lion. Mar Biol. 2009; 156: 1375–1387.

94. Pozzi L, Hodgson JA, Burrell AS, Sterner KN, Raaum RL, Disotell TR. Primate phylogenetic relationships and divergence dates inferred from complete mitochondrial genomes. Mol Phylogenet Evol. 2014; 75: 165–183. pmid:24583291

95. Cabrera A, Yepes J. Mamiferos Sud Americanos. 2nd ed. Buenos Aires: Historia Natural Ediar; 1960.

96. Harding LE, Smith FA. Mustela or Vison? Evidence for the taxonomic status of the American mink and a distinct biogeographic radiation of American weasels. Mol Phylogenet Evol. 2009; 52: 632–642. pmid:19501660

97. Lindblad-Toh K., Wade CM, Mikkelsen TS, Karlsson EK, Jaffe DB, Kamal M, et al. Genome sequence, comparative analysis and haplotype structure of the domestic dog. Nature. 2005; 438: 803–819. pmid:16341006

98. Koepfli KP, Pollinger J, Godinho R, Robinson J, Lea A, Hendricks S, et al. Genome-wide evidence reveals that African and Eurasian Golden Jackals are distinct species. Curr Biol. 2015; 25: 2158–2165. pmid:26234211

99. Viranta S, Atickem A, Werdelin L, Stenseth NC. Rediscovering a forgotten canid species. BMC Zool. 2017; 2: 6.

100. Gopalakrishnan S, Sinding MS, Ramos-Madrigal J, Niemann J, Samaniego Castruita JA, Vieira FG, et al. Interspecific gene flow shaped the evolution of the genus Canis. Curr Biol. 2019; 29: 4152. pmid:31794747

101. Higdon JW, Bininda-Emonds OR, Beck RM, Ferguson SH. Phylogeny and divergence of the pinnipeds (Carnivora: Mammalia) assessed using a multigene dataset. BMC Evol Biol. 2007; 7: 216. pmid:17996107

102. Arnason U, Gullberg A, Janke A, Kullberg M, Lehman N, Petrov EA, et al. Pinniped phylogeny and a new hypothesis for their origin and dispersal. Mol Phylogenet Evol. 2006; 41: 345–354. pmid:16815048

103. Tamura K, Kumar S. Evolutionary distance estimation under heterogeneous substitution pattern among lineages. Mol Biol Evol. 2002; 19: 1727–1736. do: pmid:12270899

104. Hassanin A. Phylogeny of Arthropoda inferred from mitochondrial sequences: strategies for limiting the misleading effects of multiple changes in pattern and rates of substitution. Mol Phylogenet Evol. 2006; 38: 100–116. pmid:16290034

105. Briggs AW, Stenzel U, Meyer M, Krause J, Kircher M, Pääbo S. Removal of deaminated cytosines and detection of in vivo methylation in ancient DNA. Nucleic Acids Res. 2010; 38: e87. pmid:20028723

106. Bromham L. The genome as a life-history character: why rate of molecular evolution varies between mammal species. Phil Trans R Soc B. 2011; 366: 2503–2513. pmid:21807731

107. Acquarone M, Born EW, Speakman JR. Field metabolic rates of walrus (Odobenus rosmarus) measured by the doubly labeled water method. Aquat Mamm. 2006; 32: 363–369.

108. Lowry L. Odobenus rosmarus. The IUCN Red List of Threatened Species 2016: e.T15106A45228501. https://dx.doi.org/10.2305/IUCN.UK.2016-1.RLTS.T15106A45228501.en.

109. Brown WM, George M Jr, Wilson AC. Rapid evolution of animal mitochondrial DNA. Proc Natl Acad Sci USA. 1979; 76: 1967–1971. pmid:109836

110. Philippe H, Forterre P. The rooting of the universal tree of life is not reliable. J Mol Evol. 1999; 49: 509–523. pmid:10486008

111. Xia X, Xie Z, Salemi M, Chen L, Wang Y. An index of substitution saturation and its application. Mol Phylogenet Evol. 2003; 26: 1–7. pmid:12470932

112. Finstermeier K, Zinner D, Brameier M, Meyer M, Kreuz E, Hofreiter M, et al. A mitogenomic phylogeny of living primates. PLoS One. 2013; 8: e69504. pmid:23874967

113. Hassanin A, Hugot JP, van Vuuren BJ. Comparison of mitochondrial genome sequences of pangolins (Mammalia, Pholidota). C R Biol. 2015; 338: 260–265. pmid:25746396

114. Gibb GC, Condamine FL, Kuch M, Enk J, Moraes-Barros N, Superina M, et al. Shotgun mitogenomics provides a reference phylogenetic framework and timescale for living Xenarthrans. Mol Biol Evol. 2016; 33: 621–642. pmid:26556496

115. Warnock RC, Parham JF, Joyce WG, Lyson TR, Donoghue PC. Calibration uncertainty in molecular dating analyses: there is no substitute for the prior evaluation of time priors. Proc Biol Sci B. 2015; 282: 20141013. pmid:25429012

116. Zachos J, Pagani M, Sloan L, Thomas E, Billups K. Trends, rhythms, and aberrations in global climate 65 Ma to present. Science. 2001; 292: 686–693. pmid:11326091

117. Zachos JC, Dickens GR, Zeebe RE. An early Cenozoic perspective on greenhouse warming and carbon-cycle dynamics. Nature. 2008; 451: 279–283. pmid:18202643

118. Kirtland Turner S, Sexton P, Charles C., Norris RD. Persistence of carbon release events through the peak of early Eocene global warmth. Nature Geosci. 2014; 7: 748–751.

119. Hooker JJ. A two-phase mammalian dispersal event across the Paleocene–Eocene transition. Newsl Stratigr. 2015; 48: 201–220.

120. Sluijs A, Schouten S, Pagani M, Woltering M, Brinkhuis H, Sinninghe Damsté JS, et al. Subtropical Arctic Ocean temperatures during the Palaeocene/Eocene thermal maximum. Nature. 2006; 441: 610–613. pmid:16752441

121. Pross J, Contreras L, Bijl PK, Greenwood DR, Bohaty SM, Schouten S, et al. Persistent near-tropical warmth on the Antarctic continent during the early Eocene epoch. Nature. 2012; 488: 73–77. pmid:22859204

122. Bininda-Emonds OR, Cardillo M, Jones KE, MacPhee RD, Beck RM, Grenyer R, et al. The delayed rise of present-day mammals. Nature. 2007; 446: 507–512. pmid:17392779

123. McKenna MC, Bell SK. Classification of mammals: above the species level. New York: Columbia University Press; 1997: 640 pp.

124. Hunt RM Jr. Basicranial anatomy of the living linsangs Prionodon and Poiana (Mammalia, Carnivora, Viverridae), with comments on the early evolution of Aeluroid Carnivorans. Am Mus Novit. 2001; 3330: 1–24.

125. Gaubert P, Veron G. Exhaustive sample set among Viverridae reveals the sister-group of felids: the linsangs as a case of extreme morphological convergence within Feliformia. Proc R Soc B. 2003; 270: 2523–2530. pmid:14667345

126. Hunt RM Jr., Tedford RH. Phylogenetic relationships within the aeluroid Carniovra and implications of their temporal and geographic distribution. In: Szalay FS, Novacek JJ, McKenna MC, editors. Mammal phylogeny. New York: Springer; 1993: pp. 53–73.

127. Veron G. Phylogeny of the Viverridae and ’viverrid-like feliforms. In: Goswami A, Friscia A., editors. Carnivoran evolution: new views on phylogeny, form and function. Cambridge: Cambridge University Press; 2010: pp. 64–91.

128. Russell DE, Tobien H. Mammalian evidence concerning the Eocene-Oligocene transition in Europe, North America and Asia. Developments in Palaeontology and Stratigraphy. 1986; 9: 299–307.