Comprehensive species set revealing the phylogeny and biogeography of Feliformia (Mammalia, Carnivora) based on mitochondrial DNA

Extant Feliformia species are one of the most diverse radiations of Carnivora (~123 species). Despite substantial recent interest in their conservation, diversification, and systematic study, no previous phylogeny contains a comprehensive species set, and no biogeography of this group is available. Here, we present a phylogenetic estimate for Feliformia with a comprehensive species set and establish a historical biogeography based on mitochondrial DNA. Both the Bayesian and maximum likelihood phylogeny for Feliformia are elucidated in our analyses and are strongly consistent with many groups recognized in previous studies. The mitochondrial phylogenetic relationships of Felidae were for the first time successfully reconstructed in our analyses with strong supported. When divergence times and dispersal/vicariance histories were compared with historical sea level changes, four dispersal and six vicariance events were identified. These vicariance events were closely related with global sea level changes. The transgression of sea into the lowland plains between Eurasia and Africa may have caused the vicariance in these regions. A fall in the sea level during late Miocene to Pliocene produced the Bering strait land bridge, which assisted the migration of American Feliformia ancestors from Asia to North America. In contrast with the ‘sweepstakes hypothesis’, our results suggest that the climate cooling during 30–27 Ma assisted Feliformia migration from the African mainland to Madagascar by creating a short-lived ice bridge across the Mozambique Channel. Lineages-through-time plots revealed a large increase in lineages since the Mid-Miocene. During the Mid-Miocene Climatic Optimum, the ecosystems and population of Feliformia rapidly expanded. Subsequent climate cooling catalyzed immigration, speciation, and the extinction of Feliformia.


Introduction
Feliformia is a large suborder of Carnivora within the eutherian clade. The monophyletic origin of Feliformia is supported by molecular and morphological data [1][2][3][4]. According to PLOS  Species 2000 and the databases ITIS Catalogue of Life [5] and Mammal Species of the World [6], Feliformia has six families (Felidae, Viverridae, Eupleridae, Nandiniidae, Herpestidae and Hyaenidae) with eight subfamilies (two in family Felidae, four in family Viverridae and two in family Eupleridae, while the families Herpestidae, Nandiniidae and Hyaenidae have no subfamilies), 54 genera and approximately 123 species. Due to pressures such as human hunting, habitat loss, lack of prey, and climate change, much of the extant Feliformia diversity is currently under extreme threat. According to The IUCN Red List of Threatened Species, Version 2016-2 [7], there are 12 endangered species: Viverra megaspila, Prionailurus planiceps, Panthera uncia, Panthera tigris, Mungotictis decemlineata, Lynx pardinus, Leopardus jacobita, Galidictis grandidieri, Eupleres major, Cynogale bennettii, Chrotogale owstoni, and Catopuma badia, and one species, Cryptoprocta spelea, has been classified as extinct. Mitochondrial (mt) DNA is a useful marker system in phylogenetic analyses because of its maternal mode of inheritance and relative lack of recombination [8]. The mt genome holds a shorter expected coalescence time compared with nuclear loci; thus, there is a greater probability that the mt gene tree will accurately reflect the species tree [9]. Although phylogenetic analyses of mtDNA have proven problematic in some cases, such as the position of snakes among squamates [10] or the relationships of rodents among mammals [11], analyses of mtDNA are often congruent with those derived from nuclear genes when appropriate sampling of taxa and phylogenetic analyses are used [12,13]. In addition, comparatively few nuclear genes are widely used in phylogenetic analyses of Feliformia; for example, only Johnson et al. [14] used multiple nuclear genes to analyze the phylogenetic relationships of Felidae. Most other Feliformia taxa lack these genes. Thus, moderately sized mitochondrial genomes remain an attractive data resource for Feliformia phylogenetics.
A phylogenetic framework is critical to discover and conserve the diversity of extant Feliformia diversity. Previous taxonomic genetic sequencing and phylogenetic analysis have revealed many subgroups of Feliformia [3,[14][15][16][17][18][19][20][21][22][23][24]. Several studies have suggested that increasing the number of sampled taxa enhances the accuracy of phylogenetic analyses (e.g., Rannala et al. [25]; Zwickl and Hillis [26]). However, Agnarsson et al. [24] evaluated the phylogenetic relationships of a comprehensive group of species of dogs, cats, and kin with only complete or partial mitochondrial Cytb gene sequences and obtained low support in many nodes of the Bayesian phylogeny tree, especially the nodes between the seven families within Feliformia. Thus, a large-scale and exhaustive species set is lacking for the phylogeny of extant Feliformia with multiple genes. A detailed historical biogeographic analysis of Feliformia also requires a comprehensive species data set to explore.
Here, we present an estimate of a comprehensive data set for the phylogeny and biogeography of Feliformia that includes 102 species (83% of the 123 known extant Feliformia species) from 52 of the 56 currently recognized genera (93%), as well as representatives from every currently delimited, extant family and subfamily base on the data matrix including up to 15,401 bp for each species from 37 complete mitochondrial genes. We also compared the divergence time with the lineages-through-time (LTT) plot and reconstructed the geographic range evolution of Feliformia with the global sea level and climate change information to discuss how global sea level change affected the evolution of Feliformia and how Feliformia responded to Cenozoic climatic change.

Taxonomic reference
Our initial taxonomy was based on the 28 September 2015 update of the Species 2000 and ITIS Catalogue of Life database [5] and was similar to that of the Mammal Species of the World database [6] at the species and generic levels, with 6 families and 54 genera. The linsang lineage (Prionodon), which recent studies have suggested is closely related to Felidae [3,14,22], was still placed within Viverridae by Species 2000, the ITIS Catalogue of Life and Mammal Species of the World [5,6]. We agreed with the classification of Prionodon by Gaubert and Veron [22] and included it in the monogeneric family Prionodontidae in our taxonomy.

Molecular data
In this study, we collected mitochondrial concatenated alignments consisting of 15,401 bp in 103 in-group species from GenBank and previous studies. These species included 83% of the currently recognized Feliformia species, representing 37 Felidae species in 15 genera (88% of 42 known extant species and 100% of 15 genera), 29 Viverridae species in 14 genera (82% of 35 species and 93% of 15 genera), 24 Herpestidae species in 12 genera (73% of 33 species and 86% of 14 genera), 6 Eupleridae species in 6 genera (75% of 8 species and 86% of 7 genera), and 4 Hyaenidae species in 3 genera (100% of 4 species and 100% of 3 genera), 2 belonging to the monogeneric family Prionodontidae and one (the only) species of the monogeneric family Nandiniidae. The concatenated alignment sequence matrix contained data from 103 species for cytb (100%, 402-1139 bp), 82 species for ND2 (79%, 696-1044 bp), and 45 species for the mitochondrial genome (exclude D-loop region) (56%, 15,401 bp). GenBank numbers are listed in S1 Table. In this study, many taxa had large amounts of missing data (some >93%), and on average each species had 50% missing cells. However, two lines of evidence suggest that these missing data are not problematic. First, two genes (Cytb/ND2) were shared by the vast majority of taxa (100% and 76%, respectively), providing a "backbone" for the placement of most taxa based on overlapping sequence data. Our phylogeny results also suggest that this sampling design can be critically important for allowing the accurate placement of taxa with extensive missing data, as opposed to having nearly the same number species with only the Cytb gene sequence [24]. Second, several recent empirical studies have shown that the supermatrix approach (with extensive missing data in some taxa) yields generally well-supported large-scale trees that are in general highly congruent with both existing taxonomies and previous phylogenetic estimates (e.g., Driskell et al. [27]; McMahon and Sanderson [28]; Pyron et al. [29]; Thomson and Shaffer [30]; Wiens et al. [31]; Agnarsson et al. [24]).

Sequence alignment, partition strategy and substitution model selection
Nucleotide sequences were aligned using Clustal X version 1.81 [32] with default parameters and manually optimized in MEGA version 5.0 [33]. Ambiguous alignments were removed using Gblocks version 0.91b [34] with the 'with half' option and default block parameters. We performed nucleotide diversity analyses based on the Cytb gene using MEGA software.
For the concatenated data set, 63 data blocks were defined (24 data blocks for RNA genes and 39 data blocks for protein genes). Each RNA gene represented a single block, and each protein gene was separated into three blocks based on codon positions. We defined 'greedy' partitioning strategies to compare partitioning strategies and selections of corresponding nucleotide substitution models under the Bayesian information criterion (BIC) using Parti-tionFinder version 1.1.1 [35]. The seven-partition scheme was chosen as the best-fitting partitioning strategy (S2 Table), and all seven partitions favored the GTR +I+G model.

Phylogenetic analyses
Three Caniformia species (Canis lupus, Vulpes lagopus and Martes flavigula) were retrieved from GenBank and used as outgroups in phylogenetic analyses (for detailed information, see S1 Table). The concatenated data set was analyzed by maximum likelihood (ML) and Bayesian inference (BI) methods separately under the seven-partition scheme. Partitioned ML analyses were implemented using RAxML version 7.2.8 [36] with the GTR +I+G model assigned to each partition. A search that combined 200 separate ML searches was performed to identify the optimal tree, and branch support for each node was evaluated with 200 standard bootstrapping replicates (-f d -b 200 option) implemented in RAxML. The partitioned BI was conducted using MrBayes 3.2 [37] under the GTR +I+G model. All model parameters were unlinked. Four Markov chain Monte Carlo (MCMC) runs were performed with one cold chain and three heated chains (temperature set to 0.1) for 10 million generations and sampled every 1,000 generations. The first 25% of the generations were discarded, and topologies and posterior probabilities were estimated from the remaining generations.

Estimates of divergence time
To further confirm the divergence time within Feliformia, another nine complete mitochondrial genomes of the following vertebrate species were retrieved from GenBank (All GenBank numbers see S1  [38] software with a relaxed uncorrelated log-normal clock model [39] was used to date Feliformia. A Yule process was employed for each of the seven partitions, and fourteen log-normal prior constraints allowing 'hard' minimum and 'soft' maximum constraints were imposed on the tree to establish the divergence time (Table 1). For the ten in-group constraints with only minimum bounds, we set the maximum bound at 55 Ma based on arithmetic medians of Feliformia and Caniformia divergence [40]. The means and standard deviations of the log-normal distribution for each calibration point were selected to ensure 95% of the probability was within the minimum and maximum limits, and the means represented the arithmetic medians of the intervals. Analyses were performed under the same partitioned substitution model with BI phylogeny with 50 million generations while sampling every 1,000th tree, and the first 10% of the sampled trees were treated as burn-in. Burn-in and convergence of the chains were determined using the program Tracer v1.6 [41].
The LTT plot was constructed using the program Tracer v1.6 and the BEAST divergence time output files. Further, we compared interesting internal nodal ages (nodal age among Feliformia families, the crown age of each Feliformia family, obvious dispersal and vicariance nodal age) and the LTT plot with historical global sea level changes reported by Haq et al. [56] and historical global climate changes reported by Zachos et al. [57] to evaluate the vicariance and dispersal history of Feliformia and its response to historical global sea level changes and climate changes.

Area delimitation and biogeographical reconstruction
We compiled the distribution data of Feliformia species from Species 2000 and the databases ITIS Catalogue of Life and the Mammal Species of the World [5,6]. The contemporary distribution ranges of Feliformia were divided into 6 areas based on their disjunctive distribution patterns and possible biogeographical barriers. Each Feliformia species was then assigned to its associated area according to its contemporary distribution range. The six areas were as follows: A, Europe-North Asia; B, South Asia; C, Africa; D, Madagascar; E, North America; and F, South America.
The ancestral area of each family of Feliformia was reconstructed using the dispersal-extinction-cladogenesis (DEC) model in Lagrange version 20130526 [58,59]. The DEC model specifies instantaneous transition rates among discrete distribution ranges along phylogenetic branches and uses the rates to access range inheritance scenarios at cladogenesis events [59]. The analysis was conducted based on the coded distribution range as defined above, and the chronogram of Feliformia species was estimated in BEAST. The maximum number of ancestral areas was limited to two, assuming that the dispersal of ancestors was similar to that of their extant descendants [60].

Phylogenetic analyses
ML and BI analyses of the concatenated mitochondrial data set produced similar topologies. the ML tree is shown in Fig 1 and  and 92% of nodes with BI posterior probability (PP) greater than 0.9. Only five nodes exhibited weak support by ML (BP <0.6) and BI (PP<0.9) or yielded different topologies in the two analyses. Our analyses support a monophyletic origin of the families Nandiniidae, Prionodontidae, Herpestidae, Eupleridae, Hyaenidae, Viverridae and Felidae and also strongly supported the family-level relationships within Feliformia with a lineage of (Nandiniidae+ ((Felidae + Prionodontidae) + (Viverridae + (Hyaenidae + (Herpestidae + Eupleridae))))), consistent with most other studies of extant Carnivora and Feliformia phylogeny [1][2][3]21]. A majority of the phylogenetic relationships within each of these families remain unchanged compared to recent estimates [3,15,19]. However, we observed significant deviations from previous phylogenies and taxonomies, which we describe below, along with proposed solutions. Within family Felidae, our results agree with other recent studies, indicating that the two monophyletic subfamilies Felinae and Pantherinae composed the family [14,61]. Six species of subfamily Pantherinae were represented in a progressive evolutionary model and were strongly supported in our analyses. Genus Neofelis and genus Panthera displayed a sister relationship, and Uncia uncia nested within Panthera, which is consistent with the results of Zhang and Zhang [61] and Johnson et al. [14]. Our analysis confirmed the phylogenetic relationship of Panthera species as (P. tigris + (Uncia uncia + (P. onca+(P. leo + P. pardus)))). Within family Viverridae, our results are consistent with other recent studies supporting the corresponding relationship at the subfamily and generic levels [3,20,22]. The four monophyletic subfamilies displayed a strongly supported lineage ((Viverrin + Genettinae) + (Hemigalinae + Paradoxurinae)). The relationships within each subfamily were as follows: Viverrin, (Viverricula+(Civettictis+Viverra));Genettinae, (Poiana+Genetta); Paradoxurinae, (Arctogali-dia+(Arctictis+(Paguma+Paradoxu rus))); Hemigalinae, (Macrogalidia+(Cynogale+(Chrotogale +(Diplogale+Hemigalus)))). Furthermore, although Prionodon was classified within Viverridae by Wilson and Reeder [6] and Roskov et al. [5], our phylogenetic trees were consistent with Gaubert and Veron [22] in redefining Viverridae as excluding Prionodon and placing Prionodon in the monogeneric family Prionodontidae. The P-distances between Prionodon and each of the other families of Feliformia were 0.177-0.191 (S3 Table), larger than the average family P-distance (0.169) of the other 6 traditional families and consistent with the classification of Prionodon as the monogeneric family Prionodontidae.
Within the family Hyaenidae, Our ML and BI analyses were consistent with Koepfli et al. [21] in that Crocuta and Proteles first obtained a sister relationship (BP = 95, PP = 0.99) and then together constituted a monophyletic family with Hyaena (BP = 100, PP = 1.0).
Within the family Eupleridae, our phylogenetic trees supported the monophyly of subfamily Galidiinae as a lineage (((Mungotictis+Salanoia)+Galidictis)+Galidia) with the crown node, and all internal nodes had 100% support according to both ML and BI analyses. However, Cryptoprocta ferox and Fossa fossana of the subfamily Euplerinae did not fit into one lineage. Three clades (Galidiinae, Cryptoprocta and Fossa) constitute the Eupleridae according to our phylogeny analyses and display a relationship of ((Galidiinae+ Fossa)+Cryptoprocta) with strong support.

Divergence times and ancestral areas
The topology of the chronogram shown in Fig 2 is nearly identical to that of the phylogenetic trees. The mean and 95% highest posterior density (HPD) intervals of our fourteen constrained nodes were all within our constraints. The crown age of suborder Feliformia is circa 43 Ma. The age of the internal nodes among families of Feliformia can be traced to the late Eocene and the middle Oligocene (Fig 2 and Table 1). The definitions of the numbers of the nodes and abbreviations of area delimitations (A-F) within Table1 are provided in Fig 2. The detail divergence information of internal nodes at the family levels, the ancestor nodes of each family and other obvious dispersal and vicariance nodes are shown in Table 1.
When comparing our new divergence time estimates with the fossil records, all means and 95% HPD intervals of divergence times at the constraint nodes obtained by BEAST are within the constraint intervals we set. Divergence times estimated among the seven Feliformia families in the present study are close to those estimated with multiple nuclear gene sequences (Eizirik et al. [62]) but considerably younger than those from one mitochondrial Cytb gene and two nuclear genes (Gaubert and Cordeiro-Estrela [3]) (Fig 3). Using multiple nuclear loci has become the standard method of resolving species relationships in difficult biological scenarios [63], and more nuclear loci used during the phylogenetic analyses lead to the phylogenetic relationship being closer to saturation [64]. Our divergence time estimated results are close to those of Eizirik et al. [62], indicating that the mt-genome is appropriate for Feliformia phylogeny analyses.
The ancestral range of subdivision/inheritance scenarios for each node of interest was initially analyzed by the Lagrange program (Table 1). Four dispersal and seven vicariance events (RP of relative nodes >0.6) are shown in Fig 2 (indicated by red branch and arrows) and were inferred as follows: (1) a dispersal event occurred from Africa (C) to Madagascar (D) between nodes 4 and 5, followed by a vicariance event between the two areas on node 5; (2) a dispersal event occurred from Africa to South Asia (B) between nodes 6 and 7, followed by a vicariance event between the two areas on node 8; (3) a vicariance event occurred between Africa and South Asia on node 13; (4) a dispersal event occurred from South Asia to Europe-North Asia (A) between node 14 and 15; (5) three vicariance events occurred between Europe-North Asia and North America (E) on node 18, 19 and 24, respectively; (6) a dispersal event occurred from Europe-North Asia to North America between nodes 23 and 24, followed by a vicariance event between the two areas on node 24. According to the Lagrange analysis, the ancestral ranges of the 7 families were determined as follows: Nandiniidae, Africa; Felidae, South Asia; Prionodontidae, South Asia; Herpestidae, Africa; Eupleridae, Madagascar; Hyaenidae, Africa; Viverridae, Eurasia.
Although some recent studies have estimated different sea level change histories using various methods (e.g., Haq et al. [56]; Haq and AI-Quahtani [65]; Miller et al. [66]; Miller et al. [67]), many of these sea level change results have the same variation trends as our dispersal event periods and vicariance event times. Thus, we only show the global sea level change history by Haq et al. [56] as representative in Fig 2. The timing of the vicariance nodes 5, 8, 13, 18, 19 and 24 was closely related to historical global sea level changes (Fig 2). The mean times of vicariance nodes 8 and 13 are located at the relative crest value during Middle Miocene sea level changes, respectively. The mean times of vicariance nodes 18, 19 and 24 are located at or near the relative trough during Late Miocene to Pliocene sea level changes, respectively. The 95% PHD intervals of these vicariance nodes also crossed at least one sea level rise or fall history. All four dispersal events have at least one rapid sea level fall during the dispersal branch times. The vicariance event between African mainland and Madagascar on node 5 was at a relatively low sea level period in the Middle Oligocene, 27.4 (24.39-31.22) Ma.

Lineages through time
The expected LTT plot of Feliformia species was constructed by the Tracer program using the output tree files from the BEAST divergence time analyses. Expected and observed LTT plots hold a similar increasing tendency, and the lineage number of the expected LTT plot was higher than the observed LTT plot (Fig 4) because the expected LTT plots were calculated under the birth-death process [68] and because there are predicted extinction lineages in the LTT plot. The number of lineages slowly increased before the middle Miocene and sharply  thereafter. Comparison of the LTT plots with historical changes in global temperature [57] revealed a strong linear increase in the number of Feliformia lineages with gradual cooling of the climate after the Mid-Miocene Climatic Optimum.

Phylogeny
Several studies have suggested that increasing the number of sampled taxa enhances the accuracy of phylogenetic analyses (e.g., Rannala et al., [25]; Zwickl and Hillis, [26]); thus, increased taxonomic sampling was performed to obtain the present data matrix. We included 86% of genera within Feliformia to successfully reconstruct the phylogenetic trees. The relationships between and within 7 monophyletic families were strongly supported by the phylogenetic analyses. Our analysis of three data subsets (Cytb gene, Cytb and ND2 genes, complete mtgenome) indicates a progressive increase in the average bootstrap support value for the six nodes among the 7 monophyletic families of Feliformia when an increasing number of genes are analyzed (Fig 5) and the complete mitochondrial genome subset, for which each Feliformia family has at least one species holding the complete mitochondrial genome sequence, shows the greatest performance according to the phylogenetic analyses.
For Herpestidae, Hyaenidae and Viverridae, phylogenetic relationships within each family were strongly supported in our analyses, and the results were consistent with previous studies [3,19,20,22]. For taxonomies within the genus Herpestes, our analyses indicated that H. ichneumon should be included within the genus Galerella and that H. naso is a close relative of Atilax paludinosus.
The relationships between all species of Euplerinae and some species of Galidiinae have been analyzed by Poux et al. [69] based on nuclear genes. Poux et al. also confirmed three major divisions within Eupleridae, including the Galidiinae clade, Cryptoprocta ferox clade, and the Fossa and Eupleres clades. Our analyses revealed similar nucleotide diversity between the three major clades (Galidiinae, Cryptoprocta ferox and Fossa fossana) with P-distances of 0.143 for Galidiinae vs. C. ferox, 0.130 for Galidiinae vs. F. fossana, and 0.142 for C. ferox vs. F. fossana. These findings are consistent with the three major divisions within the family Eupleridae.
Prionodon exhibited the same phylogenetic relationship observed in previous analyses of Felidae's closest relative but did not belong to Viverridae [3,14,22]. Our phylogeny and P-distance were consistent with Gaubert and Veron's [22] classification of Prionodon as the monogeneric family Prionodontidae.
Our mitochondrial phylogenetic analyses successfully reconstruct the phylogenetic tree for Felidae. The phylogenetic relationships between eight Felidae lineages were well resolved with the exception of the relationships between the Lynx lineage, Puma Lineage and Bay Cat lineage. The subfamily Felinae consisted of two main clades with four subclades. Our analyses clearly showed that the phylogenetic relationships between subfamily Pantherinae species were strongly represented by (Panthera tigris + (Uncia uncia + (P. onca+(P. leo + P. pardus)))). The taxonomies of Uncia uncia included a part of Panthera, which were also consistent with many previous mitochondrial phylogeny studies [14,24,61], and we agree with Yu et al. [70] that Unica is congeneric with Panthera.

Historical biogeography
Our ancestral area reconstructions show that Feliformia were distributed in South Asia and Africa before the Middle Miocene and then dispersed to America several times since the Late Miocene. Three lowland plains (the Mesopotamian Plain, the Nile Valley and the Indo-Gangetic Plain) border the three main Feliformia dispersal realms: Europe-North Asia, South Asia, and Africa. Since the Eocene, global sea levels have changed greatly and repeatedly risen by more than 100 meters [56]. The separation of these dispersal realms by drastically rising sea levels to cover the lowland plains may have contributed to vicariances and speciation events during Feliformia evolutionary history. For nodes 8 and 13, two vicariance events occurred between South Asia and Africa. The mean ages of nodes 8 and 13 are at a sea level increase period, and the 95% PHD intervals of these vicariance nodes also crossed at least one sea level rise process. Before node 8, a dispersal event from Africa to South Asia between nodes 6 and 7 occurred during two of the lowest sea levels, which were approximately 50 m higher than present, at 17-15 Ma throughout the Late Oligocene to Middle Miocene. The Eurasian and African Feliformia could spread through the lowland plains during sea level falls and separate from each other by marine transgression into the lowland plains when sea levels rise. Thus, the most reasonable vicariances of these nodes were a separation between South Asia and Africa by the transgression of the lowland plain during a period of sea level rise. Nodes 18, 19 and 24 display three vicariance events between North America and Europe-North Asia, which were separated by the Bering Strait. Before node 24, there was a rapid sea level fall during the dispersal event between nodes 23-24, followed by a vicariance with a sea level rise. The sea level fall could have facilitated neritic region exposure to provide Feliformia with the opportunity to disperse across the Bering Strait. Thus, we think that the Bering Strait provided a land bridge during rapid decreases in sea level to facilitate the dispersion of American Feliformia ancestors from North Asia to North America, followed by a vicariance between the two areas when the sea level rose.
Although the tectonic plate movement of Madagascar relative to Africa began in the middle Jurassic and ended in the Early Cretaceous [71], several studies have determined that the terrestrial vertebrates on the island of Madagascar, including amphibians and mammals, originated on mainland Africa and migrated to Madagascar between the middle Paleocene and Oligocene [69,72,73]. No evidence supports a terrestrial bridge across the Mozambique Channel. The 'sweepstakes hypothesis,' which states that ancestors of Madagascar's present-day mammalian stock rafted there from Africa, has been recognized as a valid colonization method [72,[74][75][76]. Our divergence time indicated that mainland African Herpestidae ancestors were separated from their Madagascar-based Eupleridae ancestors at circa 27 Ma (95% HPD: 24.39-31. 22 Ma) at a relatively lower sea level period that occurred during Oligocene global sea level changes (Fig 2). If the 'sweepstakes hypothesis' is valid, the relatively low sea level could have reduced the distance from mainland Africa to the island of Madagascar, thereby improving the success rate of Eupleridae ancestor over-water dispersal from mainland African to Madagascar. Moreover, we contrasted the dispersal time  of the Madagascar Eupleridae ancestor with global climate change history and found that the global climate was at an extraordinarily low temperature during this time, as estimated by many benthic foraminiferal oxygen isotopes (δ 18 O) studies [57,[77][78][79]. The continental ice sheets, perhaps even sea ice, were present at such a low global temperature climate [80]. Thus, there was likely a shortlived bridge of sea ice across the Mozambique Channel that connected Madagascar and Africa; however, we did not find any specific records or references to this. If the Madagascar-Africa sea ice existed, the migration of the Madagascar Feliformia ancestor could shorten the migration time from 25-30 days according to the 'sweepstakes hypothesis' [81][82][83] to several days, making the migration more feasible. Our results suggest that the cooling climate during 30-27 Ma may have created a short-lived ice bridge connecting Madagascar and the Africa mainland, which assisted the migration of Madagascar Feliformia out of the African mainland.

The response of Feliformia to Cenozoic climatic change
Both rapid and gradual climatic changes directly influenced the evolution and ecology of mammalian species and communities throughout the Cenozoic [84]. One of the most intriguing climatic events of the Neogene was the global warming event in the Middle Miocene (~15-18 Ma) known as the Mid-Miocene Climatic Optimum [85], after which Earth's climate gradually cooled while ice volumes increased globally [86,87]. Climatic changes during the Middle Miocene were associated with a sharp increase in the evolutionary turnover of both terrestrial and marine biota [85]. Our LTT plots also suggest that the climate since the Middle Miocene has been conducive to Feliformia diversification and speciation. The significant increase in lineages began at the Mid-Miocene climatic optimum and continued until the Pliocene. Similar increases also occurred in groupers [88]. Our hypothesis suggests that ecosystem and population expansion of Feliformia occurred rapidly during the Mid-Miocene Climatic Optimum. Subsequent climate cooling until the Pliocene, which could have directly catalyzed immigration, speciation, and the extinction of Feliformia, was a significant factor contributing to the Neogene evolutionary history of Feliformia species.