Previous studies of Mammalia phylogeny

Most studies of mammalian evolutionary history have focused on backbone-level divergences or species-level subclade radiations, but not both. For example, Carnivora (approximately 300 living species of cats, dogs, and allies; [23,51,52]) and Cetacea (approximately 90 species of whales and dolphins; [53–56]) are particularly scrutinized because of their well-studied fossils and diverse ecological habits. At the backbone level of mammalian superordinal divergences, greater paleo- to neontological integration [57–59] has recently helped bring the “rocks and clocks” of fossil-calibrated molecular ages into greater harmony [60–65]. However, controversy persists regarding both backbone node ages (e.g., [66–72]) and topological relationships (e.g., [73–75]) despite the broad application of phylogenomic and phenotypic data. Some nodes may in fact remain obstinate (e.g., due to guanine-cytosine [GC]-biased gene conversion [76,77]). Therefore, Bayesian strategies that seek to accommodate the confidence (or lack thereof) in estimated node ages and relationships, rather than collapse it to one “best” consensus, appear most valuable for testing hypotheses related to diversification processes in mammals [33,34,59,78–80].

Only a handful of studies have ever attempted to unite species-level molecular divergences with fossil ages on a Mammalia-wide basis (Table 1). The landmark study of Bininda-Emonds and colleagues [81] used a supertree approach (“matrix representation parsimony” [MRP] [38,82]) for combining source trees estimated from either DNA or morphology into a time-scaled phylogeny of 4,510 mammal species. The MRP supertree was based on the taxonomy of Mammal Species of the World, second edition (MSW2) [83] and was updated twice: (1) Fritz and colleagues [84] linked the taxonomy to 5,020 of the 5,415 species in Mammal Species of the World, third edition (MSW3) [85] and fixed errors in the dating of bats [25,86]; and (2) Kuhn and colleagues [87] resolved the >50% of unresolved nodes (2,503 polytomies) remaining in the MRP supertree using a stochastic birth–death model, creating a set of 1,000 trees with random variation in the placement of unresolved branches [87]. Versions of the MRP supertree have been widely applied to questions of species diversification (e.g., [1,9,88,89]) and conservation (e.g., [25,29,90,91]) despite the initially unresolved species and consequent potential for artifacts in downstream analyses, in part because it contained the only estimates of evolutionary branch lengths across most of Mammalia.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Previous species-level phylogenies of Mammalia relative to the present study.

https://doi.org/10.1371/journal.pbio.3000494.t001

Only two other estimates of species-level mammalian phylogeny have been published. One was the DNA-based supertree of Faurby and Svenning [92], and the other was the consensus timetree of Hedges and colleagues [93] (Table 1). The DNA supertree was constructed from the hierarchical merging of 290 overlapping subtrees (mostly estimated at the level of genus or family), followed by random resolution of polytomies, addition of DNA-lacking species (n = 2,364), and subsequent rescaling to time using mean node ages from secondary (e.g., [60]) or tertiary (e.g., [81]) sources. This study was an improvement over the previous MRP supertree by directly estimating Bayesian subtrees using DNA sequence data (range: 1–26 markers for most species, up to 56 for Primates), which allowed greater phylogenetic uncertainty to be included in their final distribution of 1,000 trees. Similarly, the expansive “timetree of life” [93] is also a supertree, albeit on the much larger scale of eukaryotes initially and then pruned to Mammalia for application to several subsequent rate-based phylogenetic analyses (e.g., [96–101]). Hedges and colleagues [93] merged the overlapping topologies and mean node ages from 91 divergence-time studies of mammals published from 1991–2013 (see http://www.timetree.org/references), using an approach they call “hierarchical average linkage” to construct the supertree topology. The Hedges and colleagues [93] timetree is a single consensus estimate (one tree). Because supertree algorithms produce polytomies when overlapping sources disagree, the corresponding branch lengths were secondarily time scaled or simulated in these studies as well as for the MRP supertree. Thus, evolutionary rate estimates in these supertrees are expected to be unreliable, as the authors of the DNA supertree admit: “Our approach places the greatest weight on the topology, which means that analyses using the resulting phylogeny should focus on the topology rather than on branch lengths” (see [92], page 16).

A path forward

Reconstructing species-level mammal trees has forced researchers to depart from the standard phylogenetic approaches for jointly inferring species relationships and node ages from primary character data (molecular or morphological; reviewed in [102]). Steps of merging overlapping sources, collapsing conflicting nodes, and applying point-estimate dates to scale phylogenies to time are common to the MRP supertree, DNA supertree, and consensus timetree analyses. In each case, branch length information is reduced from the original data as a trade-off for enabling large-scale inference. However, with increasing computational ability and growing public databases of DNA sequence information, supertree methods are no longer the only way to infer big trees. We here leverage computational power from the Cyberinfrastructure for Phylogenetic Research (CIPRES) Science Gateway [103] and extensive public data deposited in the United States National Center for Biotechnology Information (NCBI) (Genbank [104]) to filter, clean, assemble, and then reconstruct phylogenetic history from an inclusive DNA supermatrix of mammalian species. As we outline below, the goal of jointly inferring tree topology and node ages is now computationally feasible for Bayesian analysis of large clades (800–1,000 species), opening the door for greater resolution of macroevolutionary tree shape in mammals and other taxa.

DNA sequence alignment, gene trees, and error checking

Our BLAST-based [105] pipeline initially yielded 209,294 matching hits across all 31 genes. We used an iterative per-gene approach to clean annotation errors in NCBI (Fig 2A), as follows: (1) sequence alignment, (2) error checking for stop codons and insufficient alignment overlap, and (3) gene tree construction (Randomized Axelerated Maximum Likelihood [RAxML] v.8.2.3 [49]—see S1 Text, section 3). To minimize stop codons for the 26 coding fragments (mitochondrial DNA [mtDNA] and exons), we aligned each to the appropriate amino acid reading frame and excluded unaligned (entirely nonoverlapping) sequences, as well as rogue taxa (see S1 Text, section 3). In total, our error-checking steps excluded 1,618 sequences across all genes (i.e., 7.2% of the 22,504 individual DNA sequences after taxonomic reconciliation; S2 Table). These exclusions corresponded to 119 species, yielding 4,098 species with ≥1 gene fragment validated in the final 31-gene matrix (S1 Data lists excluded sequences). Our procedure of DNA-baits searching, curation, and alignment of sequences from the NCBI database resulted in taxon sampling that ranged from 191 to 3,581 species per gene (Table 2).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Alignment details for the 31 gene fragments used in this study.

https://doi.org/10.1371/journal.pbio.3000494.t002

The 31-gene supermatrix

Concatenation of the per-gene alignments was performed in Geneious v.9.1 [106], resulting in a sites-by-taxon supermatrix of 39,099 base pairs (bp) and 4,098 species that was 11.9% complete in terms of ungapped sites. The final DNA supermatrix consisted of 21,021 DNA sequences from public databases (see S8 Fig for top individual contributors). We evaluated partitioning schemes for the supermatrix using PartitionFinder v.1.1.1 [107], finding that a nine-partition model was most suitable (Table 3). This model has a combined partition for APP, CREM, and FBN1 and then one partition each for BMI1; PLCB4; and first, second, and third codon partitions for nuclear DNA (nDNA) exons as well as for mtDNA fragments. For all partitions, either general time-reversible plus gamma (GTR + G) or plus gamma and invariant sites (GTR + I + G) was the best model of nucleotide evolution. We chose the simpler GTR + G model for all downstream phylogenetic analyses because including both I and G types of rate heterogeneity is known to make both model parameters difficult to estimate [49,108].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Model parameters optimized during the global RAxML tree search of the 4,098-species supermatrix.

https://doi.org/10.1371/journal.pbio.3000494.t003

Global RAxML tree

Phylogenetic analysis of the 4,098-species DNA matrix was first performed in RAxML with the goal to identify the single best-supported topology for mammals. For RAxML, we ran five independent analyses, each specifying 100 bootstrap replicates and using the “-f a” option and GTRCAT model to search for the best-scoring tree using ML (with this setting, the ML optimizations start from every fifth bootstrap tree [49]). Each RAxML analysis took approximately 5.7 days on 12 nodes of four threads each on the XSEDE cluster (Extreme Science and Engineering Discovery Environment; accessed via the CIPRES Science Gateway [103]). We subsequently summarized this single best ML tree (likelihood −3,383,607.6, tree length 255.3) by rooting it with Anolis and annotating nodes with bipartition values from 100 bootstrap replicates (S3 Data).

Patch clade delimitation

We divided the mammalian phylogeny into 28 patch clades that were nonoverlapping in their in-group species membership (Fig 2C; Table 4). Criteria for delimitation were clade size, evidence for monophyly (in our global ML tree and previous studies), and the structure of interclade phylogenetic relationships. Nodes with >75% bootstrap support (BS) were deemed well supported. The main challenge was to balance reasonable assumptions of monophyly with maximum patch size (number of species) for which we could feasibly perform Bayesian joint estimation of topology and branch lengths in less than 1 month. If Markov chain convergence were to take longer than 1 month, the need to iteratively conduct sensitivity tests and model tuning would have been unreasonable. We used MrBayes v.3.2.6 [50] for all Bayesian inference of patch clade and backbone phylogenies because of its flexible application of topological constraints. By experimentation, we concluded that approximately 800 species was the feasibility limit for patch clade size, although matrix size and complexity also influenced run times. Near this maximum, our largest patch clade (Muridae, 778 species) took 3.7 weeks to finish 33,330,000 generations in parallel on 16 BEAGLE (Broad-platform Evolutionary Analysis General Likelihood Evaluator)-enabled compute nodes. With MrBayes run times of 1.5 to 4.5 weeks for clades >200 species (Table 4), we estimate that approximately 80 weeks of run time was applied to the DNA-only and completed patch analyses for a total of 215,040 cpu hours (560 days * 24 hours * 16 nodes for final models, not counting troubleshooting).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Summary of patch clade and matrix composition.

https://doi.org/10.1371/journal.pbio.3000494.t004

Delimiting appropriate patch clades was especially challenging in bats and rodents. Here, species richness is highest, but so is missing genetic data and topological uncertainty. In the mouse-related clade of rodents (1,768 total species, 64% genetic sampling; [109]), we addressed this issue by dividing data into two large and likely monophyletic clades (Muridae and Cricetidae) and several smaller Muroidea patch clades for which interrelationships are uncertain (Dipodidae, Spalacidae, Nesomyidae, Calomyscidae, Platacanthomyidae; [110]). We thus avoided assuming a backbone topology for mouse-related rodents; instead, uncertainty in patch interrelationships was captured on the dated Mammalia backbone (see below). Note that these smaller patches were each well supported in the global ML tree except the Nesomyidae of Madagascar (BS 72), for which monophyly is well supported in other studies [110,111].

For bats, major topological uncertainty lies within Yangochiroptera (902 species, 67% genetic sampling), especially among its most basal divergences [112–114]. However, the compute time required to run Yangochiroptera as a single patch clade was prohibitive (initial attempts suggested 6–8 weeks) and with no guarantee of convergence (matrix <10% complete). Rather, we divided this group (94 BS value, S3 Data) into three patch clades:

1. Noctilionoidea (Phyllostomidae, Mormoopidae, Noctilionidae, Thyropteridae, Furipteridae, Mystacinidae, and Myzopodidae);
2. Vespertilionoidea (Vespertilionidae, Molossidae, and Natalidae); and
3. Emballonuroidea (Emballonuridae and Nycteridae).

Most controversial of these delimitations is the placement of Myzopodidae, which we include with Noctilionoidea according to BS 76% in the global ML tree (alternatively linked to Emballonuridae [113]). Support for the Vespertilionoidea was uncertain in our global ML tree (BS 51 joining Natalidae with Molossidae + Vespertilionidae; BS 48 for Vespertilionidae with Mollosidae), and similarly, we recovered Emballonuroidea with BS 52. Nevertheless, our patch clade schema represents the best-supported hypotheses for Yangochiroptera, with diversity divided into manageable group sizes.

The remaining patch clades encompassed not only major swaths of mammalian diversity (e.g., Marsupialia, Primates) but also very small clades like Monotremata (5 species), Pholidota (8), and Dermoptera (2; Table 4). The structure of the phylogeny and backbone uncertainty necessitated small clades to minimize unsupported monophyly assumptions. Our smallest patch clades (Dermoptera and Platacanthomyidae) were needed for this reason—however, because phylogeny estimation requires at least four taxa, we added two in-group species for MrBayes runs (Callithrix jacchus and Gorilla gorilla, and Rattus norvegicus and Spalax ehrenbergi, respectively). These species were pruned out before rescaling and pasting to the backbone.

We followed Jetz and colleagues [27] in classifying DNA-sampled species as type 1 (sampled for one or more genes) and DNA-missing species as type 2, 3, or 4, as follows: type 2, DNA available for at least one congeneric species (constrain to genus); type 3, no DNA in the genus, but available in the same family (constrain to family); and type 4, no DNA in the family, but available in the same order (constrain to order).

Along with 4,098 type 1 species, we had 1,649 species in the type 2 category, meaning that 91% of the 1,813 DNA-missing species could be constrained to a DNA-sampled genus. Beyond that, we had 115 genera entirely unsampled for DNA, to which 156 type 3 species belong. Most of these missing genera are rodents (73 genera, 58 of which are muroids) or bats (22 genera). Additionally, there were three extinct families in our taxonomy to which no DNA was available at the time of download (eight species in the type 4 category): Nesophontidae, Prolagidae, and Thylacinidae (S1 Data).

Tip versus node dating on the mammalian backbone

Comparing our node- and tip-dating analyses, we find broadly similar backbone ages with a few exceptions (Fig 3A). Chiefly, tip dating produced an older root of Mammalia and younger divergences among some rodent and bat lineages than did the node-dating analyses (Fig 3A; S3 Table). Tip dating posits that crown mammals began radiating as early as the Permian–Triassic boundary (Fig 2B; approximately 246 Ma [222.1, 268.3]), but this is much more likely to have been an early Jurassic event [116–120]. Other recent tip-dating analyses have also recovered old ages for the Mammalia crown (e.g., approximately 204 Ma in Lee [121]), suggesting that tip dating may require a combination of root age constraints [72] and a fossilized birth–death (FBD) prior that accounts for nonrandom (diversified) taxon sampling [64] to bring estimates closer to the strict fossil age of approximately 166 Ma [66]. Here, we used the latter but not the former.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Tempo of mammalian divergences recovered in the two dating analyses of the present study.

Tip-dated (fossilized birth–death) and node-dated (exponential priors) analyses yielded broadly similar results. (a) Among major clades, mean divergence times and 95% highest posterior density intervals are shown for the 28 backbone lineages present in the full trees. (b) Species-specific rates of speciation were compared using the tip DR metric, as calculated upon 10,000 trees as harmonic mean estimates (colored dots by higher taxon) and 95% CIs (Spearman’s r = 0.93 of tip- to node-dated harmonic means). Dryad data: https://doi.org/10.5061/dryad.tb03d03; phylogeny subsets: http://vertlife.org/phylosubsets. CI, confidence interval; tip DR, tip-level pure-birth diversification rate.

https://doi.org/10.1371/journal.pbio.3000494.g003

Although tip dating recovers an older root age, we find it yields younger ages than node dating for divergences between and among Muridae, Cricetidae, and Nesomyidae in the mouse-related clade and between Noctilionoidea and Emballonuroidea in Yangochiroptera bats (inset in Fig 3A; S4, S5 and S6 Figs for further details). These same areas are topologically uncertain in both backbones, indicating that the lack of monophyly constraints in the tip-dating analysis (versus 18 in the node-dating analysis) is influencing node ages. Hence, resolving the topology of difficult nodes in the rodent and bat radiations is a question deeply intertwined with resolving their divergence times. Greater applications of phylogenomic data (e.g., [122]) as well as methods that explicitly account for life-history biases among lineages (e.g., CoEvol [62]) are promising strategies toward those joint temporal and topological goals.

Overall, tip dating is laudable for its probabilistic placement of fossils using morphological synapomorphies relative to living taxa because doing so requires fewer “hard” assumptions of fossil crown-versus-stem placement [58,123]. However, for the reasons outlined above, we have more confidence in our node-dating analysis. The indirect use of fossil data as node priors also remains more mainstream (e.g., [69,124–128]). We thus focus discussion on how the node-dating results influence the Mammalia-wide trees relative to previous studies.

Backbone-level divergences

Comparing our study with previous fossil-calibrated molecular trees reveals a growing consensus for the tempo of superordinal divergences in mammals (Fig 4). We find broad agreement (overlapping 95% confidence limits) for the crown age of Marsupialia in our study (approximately 79 Ma; 67.9, 92.8) relative to 68–97 Ma in previous studies (Fig 4; [60,61,81]). Similarly, for Placentalia, our crown estimate of approximately 92 Ma (77.4, 105.0) is concordant with previous studies including the tip-dating study of Ronquist and colleagues ([64]; approximately 85 Ma [76, 93]—but note an older placental age with a different tip-dating tree prior [72]; approximately 132 Ma [119, 148]). The consensus interpretation of the fossil record as given by Foley and colleagues [63] gives a wide allowance for the placental crown to be at least 65.2 Ma (Purgatorius stem primate [129]) and no older than 131.5 Ma (Eomaia stem eutherian [130]; arguably to Juramaia at approximately 160 Ma [117]). Nevertheless, the strict fossil-based perspective for marsupial and placental crown ages fixed at 64.85 Ma [66] appears untenable given joint consideration of the molecular and fossil evidence.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Mammalian backbone-level divergence times in our study (node-dated analysis) versus previous studies and the fossil record.

The right-side phylogeny depicts relationships among the 27 extant orders (labeled in capital letters and nested in a hierarchical list), and the dotted line represents the K-Pg extinction event, 66 Ma. Divergence times are colored per study as mean ages and 95% confidence intervals. Fossil-calibrated molecular ages are compared with min and max ages for the oldest crown fossil according to Foley and colleagues [63] and oldest stem fossil according to the Paleobiology Database. Asterisks (*) on taxon names denote three instances of “zombie lineage” disagreement of our study with previous interpretations of the fossil record (see Discussion). Note that extant Microbiotheria and Tubulidentata are monotypic, and so they lack crown ages. Dryad data: https://doi.org/10.5061/dryad.tb03d03; phylogeny subsets: http://vertlife.org/phylosubsets. K-Pg, Cretaceous–Paleogene; Ma, million years ago; max, maximum; min, minimum. Artwork from phylopic.org and open source fonts (see S1 Text, section 9 for detailed credits).

https://doi.org/10.1371/journal.pbio.3000494.g004

Particular controversies exist regarding whether early divergences in crown placentals occurred before, after, or during the Cretaceous–Paleogene (K-Pg) mass extinction event, 66 Ma [131]. Here, we recover the first four placental divergences unambiguously preceding the K-Pg: (1) Atlantogenata (when present; see below), (2) Boreoeutheria, (3) Laurasiatheria, and (4) Euarchontoglires (Fig 4). The next 21 divergence events subsequently have confidence limits that overlap the K-Pg, including 12 superordinal divergences and nine of 18 crown orders (Fig 4; S3 Table). The K-Pg event being possibly concurrent with nine of the 18 placental orders compares with previous studies finding three [61], five [81], or six [60] orders with K-Pg-overlapping divergence times (S3 Table). Our finding that no placental crown ordinal radiation definitively preceded approximately 66 Ma counters previous evidence that Eulipotyphla [60,81] and possibly Rodentia and Primates [64,81] began radiating before this event (Fig 4).

Our node-dating results are conservative with respect to the fossil record. Considering the oldest fossil genus per extant mammalian order (yellow squares in Fig 4; data from the Paleobiology Database [132]), we found consistent agreement with the expectation for these fossils to be members of the stem lineage and thus older than the crown ages from the phylogeny. These maximum fossil ages, when available, are found to be either older or overlapping our divergence-time intervals in all but two cases (Fig 4). These exceptions are (1) Diprotodontia, in which the fossil genus Paljara (Pseudocheiridae) may be as old as approximately 34 Ma [133] versus approximately 49 Ma for the crown order in our tree; and (2) Eulipotyphla, in which Litolestes and Oncocherus are Erinaceidae from as old as approximately 62 Ma [134] versus approximately 75 Ma in our tree. In both cases, molecular ages extend back further, suggesting that those fossils are either legitimate crown rather than stem members of those orders or are later-surviving stem representatives. The former case is supported for Litolestes ([66]; S4 Table compares these fossils to our stem ages).

An additional check relative to the fossil record was to search for “zombie lineages” [125], in which molecular divergence dates are younger (more recent) than the minimum ages implied by well-supported crown fossils. Comparing our dates with the consensus node calibrations of Foley and colleagues ([63] updated from those of [60]), we find broad agreement but three notable exceptions (asterisks on taxon names in Fig 4). First, and most substantially, our crown age for Sirenia (manatees, dugongs, and sea cows; approximately 13 Ma: 7.0, 22.6) is reconstructed as younger than the minimum age constraint of 41.3 Ma given in Foley and colleagues [63] and, thus, “undead” for at least 20 Ma. We reconstruct the sirenian stem divergence at approximately 54 Ma (41.5, 67.3), which implies a long stem to the crown divergence of Dugong, Trichechus, and Hydrodamalis—rather than the perspective in which those three modern genera are deeply divergent from each other [60,63,135]. This issue hinges on the acceptance of the fossils Halitherium and Eotheroides as crown sirenians because they form the minimum age constraint in previous studies. The most recent cladistic analysis of Springer and colleagues [135] found 40% BS for the placement of Halitherium and Eotheroides inside of crown Sirenia (stem taxa of the Dugong–Hydrodamalis clade, to the exclusion of Trichechus [135]). Based on our criteria for fossil inclusion [57], these fossils were placed too tenuously for use as Sirenia crown constraints. Instead, we relied on a single node prior for Afrotheria (calibration 7 in S1 Text) and, thereby, placed greater weight on molecular evidence for this node. Philosophically, we recognize that we set a high bar for placing fossils using cladistic analyses, but we contend this approach is necessary to avoid false confidence regarding the timescale of mammalian evolution.

Second, molecular divergences within the sister clades Perissodactyla and Artiodactyla in our analyses also display apparent zombie tendencies relative to some interpretations of the fossil record (Fig 4). The crown or stem placement of constraint fossils is again in question. We recover the crown divergence of Perissodactyla at approximately 39 Ma (32.6, 45.0), which overlaps the age estimate obtained in dos Reis and colleagues [61] of 52.6 Ma (41.8, 61.0) and is similar to the mean age estimates of Phillips [126] of 41.4 Ma and 36.1 Ma based on strict and relaxed molecular clocks. In contrast, calibrating Perissodactyla with the fossil genus Hyracotherium sets an age range of 55.5–61.6 Ma [60,63], which was closely mirrored in the Perissodactyla age of 56.8 Ma (55.1, 61) recovered in Meredith and colleagues ([60]; Fig 4). However, the only two cladistic studies of Hyracotherium show it falling outside of crown Perissodactlya: O’Leary and Gatesy [136] and Spaulding and colleagues [137]. Both studies recovered Hyracotherium as stemward to the clade that includes crown Perissodactlya + Artiodactyla; therefore, this fossil is actually two nodes back from being able to serve as a crown Perissodactlya constraint. The next candidate fossil for the oldest crown Perissodactyla is younger than the fossil we used to calibrate Artiodactyla: Himalayacetus subathuensis from the early Eocene approximately 52.4 Ma of India, which is the oldest stem whale according to the cladistic analysis of O’Leary and Uhen [138] (calibration 16 in S1 Text, following the compendium of Benton and colleagues [139]).

Our use of Himalayacetus to calibrate Artiodactlya, in turn, informs our recovered age of approximately 39 Ma (32.7, 46.4) for crown Whippomorpha (whales + hippos). Foley and colleagues [63] use that same fossil as a crown constraint for the Whippomorpha node (52.5–61.6 Ma), which is three full nodes tipward from the Artiodactyla crown, where we used it. Clearly, this is another case of differently interpreting the fossil record. Himalayacetus is known only from a partial dentary and two molars [140] and is tentatively allied as a stem whale [66] but is more conservatively a stem whippomorphan for use in calibrating the Ruminantia–Whippomorpha node as Benton and colleagues [139] recommend (we did not do this to avoid using the same fossil twice). Although our crown ages for Whippomorpha and Cetacea are somewhat young, the former overlaps previous age estimates of 52.2 Ma (41.9, 62.6; [81]) and 48.1 Ma (45.9, 50.1; [61]), and the latter is congruent with three previous studies (see Fig 4).

In summary, interpretations of the fossil record that lose sight of the need to cladistically confirm the placement of calibration fossils inside the crown node they are constraining appear to cause the three putative zombie disagreements between our study and Foley and colleagues [63]. Our more conservative application of the fossil record aims to exclude opinion-based assignments of fossils in crown clades [57]. The fossil record provides critical data about past mammalian diversity [11,141,142], but because the preservation of mammal fossils is spatially, temporally, and taxonomically biased (e.g., Cenozoic of North America [143]; bats [144]), we contend that providing greater weight to the molecular data is warranted. It is possible to place a priori constraints on the ages of nearly all nodes in the mammalian backbone (e.g., 84 of the 163 nodes in Meredith and colleagues [60]). However, doing so requires the use of taxonomic opinions to place fossils relative to given crown groups. As we detailed above, some of these opinion-placed fossils are subsequently found to be stemward of the calibrated node in cladistic analyses. Because both “clocks” and “rocks” have shortfalls [59,145], the implementation of fossil calibration approaches to molecular data should aim to propagate age uncertainty rather than overly restrict it, thereby enabling conservative tests of evolutionary history and its causal underpinnings.

Backbone-level topology

The objective of our study was to provide novel resolution on the rates and timing of mammalian divergence events, but these results are nevertheless relevant to a few long-standing issues of topological relationships among major clades (see S1 Text and S4, S5 and S6 Figs for detailed comparisons of the backbone consensus trees). We highlight four regions of the placental backbone that are especially controversial:

1. The rooting of Placentalia. We recover support of 0.53 posterior probability (PP) in favor of the Atlantogenata rooting (Xenarthra + Afrotheria) compared with 0.47 PP for the Afrotheria rooting (Exafroplacentalia) in the ND analyses (Fig 2B and 2C), whereas the tip-dated backbone recovered the Afrotheria rooting most commonly (0.44 PP; rooting of Atlantogenata was also recovered). The high uncertainty we recover for this basal divergence is typical of other molecular studies [60,61,63,73,75,128], although the Atlantogenata rooting has received more support in phylogenomic data sets (e.g., [75]). In contrast, studies that filter genes based on their likelihood of incomplete lineage sorting (proxied by adenine-thymine [AT] content; [73,77]) generally favor the Afrotheria rooting.
2. The position of treeshrews (Scandentia) relative to colugos (Dermoptera) and Primates. We find treeshrews allied with colugos (0.78 and 0.84 PP in node- and tip-dated analyses), and that clade is always adjacent to Primates. By comparison, Scandentia has varied in position considerably depending on analysis methodology in other studies, mostly between the result we recovered and rooting outside all other Euarchontoglires (including rodents and lagomorphs; e.g., [73,75,128]).
3. The position of guinea pig–related rodents (Hystricomorpha, also called Ctenohystrica; see [146]) relative to mouse- and squirrel-related clades. We find this controversial node, which was formerly questioned to even be inside Rodentia [147], to be unequivocally recovered as ([guinea pig, squirrel] mouse) in all backbone analyses. Strong support for this relationship was recovered in some studies [60,75,122], but others have supported squirrels outside other rodents, either with a total-evidence approach [148] or when taxon sampling is smaller [61,73,128]. Transposon evidence suggests that ancient hybridization may be complicating the early history of rodents [109], as might the disparate rates of molecular evolution in these three clades [122]. Regardless of the order of branching, these basal rodent divergences were very rapid and possibly even simultaneous (i.e., overlapping error bars for nodes 45–47 in S6 Fig).
4. The branching order of mouse-related rodent families (backbone of Supramyomorpha [146]). We recover the infraorder Myomorphi as marginally sister to Anomaluromorphi to the exclusion of Castorimorphi (0.52 PP in the node-dating backbone), as well as Muridae–Cricetidae to the exclusion of Nesomyidae somewhat favored (0.57 PP; Fig 2B). These nodes were also equivocal in the largest phylogenomic data set yet leveled at this question [122], and the study of Steppan and Schenk ([110]; six loci, 904 muroid taxa) found 93% ML BS for the Muridae–Cricetidae relationship. Again, these interfamilial mouse-related divergences appear to have been extremely rapid.

We emphasize the retained uncertainty in the placental backbone divergence (Fig 2C) as a strength of the backbone-and-patch approach because having two levels of nonoverlapping Bayesian analysis enables temporal information to be passed forward to the species tips. Rather than selecting one “best” topology for rooting the placental radiation, our trees propagate the implications of both Atlantogenata and Afrotheria rootings (and other uncertainties) to the final species-level sets of 10,000 trees, providing investigators with a more realistic resource for hypothesis testing than any single tree alone.

Species-level tree shape

Comparing the accumulation of lineages through time demonstrates how the temporal and phylogenetic uncertainty propagated in our mammal trees surpasses that of previous species-level studies (Fig 5). Viewing these phylogenies on a Mammalia-wide basis (Fig 5A), we see that the range of ages incorporated in the original MRP supertree [81,84], and subsequent analyses that resolved polytomies [87], is considerably narrower than the ages encompassed in our credible tree sets. The large number of polytomies originally present in the MRP supertree is shown concentrated at approximately 50 Ma and approximately 30 Ma (Fig 5A), particularly within the rodents (Fig 5B, lower right, top line). Resolving those polytomies changes the tree shape but does not reflect the considerable uncertainty in node ages and relationships. That is, the unresolved nodes produced in supertree studies when nodes conflict are “soft” polytomies, in which the data needed to resolve a given node is lacking [41], as opposed to “hard” polytomies, in which historically rapid divergence has led to a star phylogeny [149]. Collapsing uncertainty into soft polytomies was a purposeful tool for supertree methods to yield a single consensus picture of evolutionary topology for more species than possible under joint inference [38,41,150].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Accumulation of lineages through time for living mammal species as reconstructed in our study versus previous studies.

(a) The shape of Mammalia-wide phylogenies is compared among studies using the natural log of lineage accumulation (see legend colors). Some studies produced one consensus tree (single line), whereas other studies produced sets of 1,000 or 10,000 trees (many lines), in which case 100 trees were randomly sampled. (b) Each of the main species-level Mammalia studies is compared for three major placental orders: Rodentia (purple), Chiroptera (red), and Primates (orange). The degrees of phylogenetic uncertainty present in the tree sets is represented by the width of the lineage accumulation curves. The gray lines in the lower-right-side plot pertain to the MRP supertree with polytomies, whereas the colored lines result from randomly resolving those polytomies into a set of 1,000 trees, of which 100 trees are plotted here. Dryad data: https://doi.org/10.5061/dryad.tb03d03; phylogeny subsets: http://vertlife.org/phylosubsets. Ma, million years ago; MCC, maximum clade credibility; MRP, matrix representation parsimony. Artwork from phylopic.org and open source fonts (see S1 Text, section 9 for detailed credits).

https://doi.org/10.1371/journal.pbio.3000494.g005

The danger, of course, has been when soft polytomies are misinterpreted by subsequent investigators who assume that all temporal and phylogenetic signatures in supertrees are driven by biological processes. For example, the study of Stadler and colleagues [88] made an important modeling advance for detecting tree-wide shifts in diversification rates, but the biological conclusion of a major rate shift approximately 30 Ma in rodents was apparently driven by soft polytomies in the MRP supertree (see Fig 5B). Miscommunication between the stated purpose of supertrees—“to produce phylogenies based on all data sources” [150: 266]—and the need for big trees to additionally model all uncertainty in those data sources appears to have limited the durability of supertree-based inferences and, perhaps, non-Bayesian methods generally [33].

The other mammal supertrees similarly contain less temporal uncertainty in their lineage accumulation curves than the backbone-and-patch trees of this study (Fig 5). Constructed directly from genus- and family-level DNA trees, the DNA supertree study of Faurby and Svenning [92] represents an advance over the MRP supertree. However, as shown in the lineage accumulation curves of rodents, bats, and primates (Fig 5B, top right), there are unusual artifacts of limited temporal uncertainty between those crown ordinal divergences until approximately 55 Ma, when the curves broaden to represent greater rate uncertainty. Are paleomammalogists actually more certain about the timing of events near the K-Pg extinction event 66 Ma than they are about modern divergences? Although this seems unlikely given preservation biases in the fossil record (e.g., [151,152]), that is the information conveyed by the DNA supertree. However, rather than an intended statement of confidence by the study’s authors, this is an artifact of the hierarchical merging and rescaling of overlapping subtrees onto a time-scaled backbone that lacks age uncertainty. Although the authors of the DNA supertree contend that applications “should focus on the topology rather than on branch lengths” (see [92], page 16), that advice may be subsequently ignored by researchers because of the allure of addressing important questions in comparative biology. Parsimony-based ancestral state reconstruction is perhaps the only methodological approach that entirely neglects branch lengths [102,153]. In this context, nearly all phylogenetic questions are to some extent “rate based,” although the relative importance of tree- and rate-based information to different questions is subject to debate. The key point is that researchers seeking to perform all but the most basic parsimony analyses are aided by phylogenies that propagate uncertainty in both rates and topology.

Tip-level speciation rates of mammals

We calculated tip-level speciation rates across all living mammal species (Fig 1) for comparison with those estimated on supertree phylogenies (Fig 6). We use the tip-level pure-birth diversification rate (tip DR) metric [27] because it is readily calculable across all 10,000 trees in our credible sets while being highly correlated with model-based estimators of tip speciation rates (demonstrated in Quintero and Jetz [154] and reviewed in Title and Rabosky [28]). The reciprocal of the tip DR metric is a statistic called “equal splits” [24], which is tightly related to the “fair proportion” statistic commonly used to determine evolutionary distinctiveness (ED) (e.g., [39,40]). However, the ability to robustly estimate tip DR and ED requires trees that are completely sampled (contain all modern species) and probabilistically inferred (with uncertainty in topology and branch lengths).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Tip-level speciation rates (tip DR metric) from this study versus previous Mammalia-wide phylogenies.

Comparisons (a) as plotted on trees with relative color scales calibrated per data set so that the top 1% of the tip rate harmonic means correspond to the bright red color of each tree and (b) on a pairwise basis for all species with taxon names matching directly between data sets (n = 4,670, 5,329, and 5,033 species, respectively). Note that the x-axes differ but correspond to the range of tip DR values (95% CI) of each data set. Dryad data: https://doi.org/10.5061/dryad.tb03d03; phylogeny subsets: http://vertlife.org/phylosubsets. CI, confidence interval; tip DR, tip-level pure-birth diversification rate.

https://doi.org/10.1371/journal.pbio.3000494.g006

Broadly, we find substantial heterogeneity in tip rates across the mammal tree, sometimes with a few high-tip-rate species nested among low-tip-rate species (Fig 1), resulting in long right-side tails in the tip rate distributions (positive skew, e.g., clades 38 and 44 in Fig 1). We find the consistently highest tip speciation rates in simian primates (clades 42–43 in Fig 1), including the human genus Homo (80th percentile, median 0.321 species/lineage/Ma; Homo sapiens and three extinct species) and Indo-Malayan lutung monkeys (95th percentile, 0.419, Trachypithecus). In contrast, species of Ctenomys tuco-tucos and Pteropus flying foxes display high tip speciation rates among otherwise slower-evolving species (clades 44 and 38, Fig 1). The evolutionarily distinctive platypus and aardvark have the lowest tip speciation rates (clades 1, 14; Fig 1). We suggest that tip rate skew is measuring aspects of within-clade speciation rate variation that may be otherwise uncaptured by model-fitting approaches (S5 Table). Future studies may thus find clade-level distributions of tip rates to be useful for comparative analysis.

Assessing how the different temporal frameworks of the node- and tip-dated backbones influence the species-level rate calculations (Fig 3B), we find that mammal species have broadly similar tip DR estimates across our tree sets. Indeed, there is approximately the same amount of variation in the 95% CIs of tip DR within a given tree set as between the two sets. The tip-dated phylogenies produce somewhat higher estimates of the tip DR harmonic mean (maximum of approximately 1.1 species/Ma versus approximately 0.8 in the ND phylogenies) but are nevertheless strongly correlated to the ND estimates (Spearman’s r = 0.93; linear model: y = 0.02 + 1.17x, R² = 0.85). The majority of the variation in tip rates among phylogenies appears to trace back to the younger node ages for mouse-related rodents and yangochiropteran bats (Fig 3A). Nevertheless, the internal consistency of each tree set suggests that applying either the tip- or ND phylogenies (or both) to a given comparative analysis would be appropriate.

We next compared the tip DR values from our backbone-and-patch analysis to estimates for the same species in previous supertrees to understand how different tree-building methods influence those rates (Fig 6; tree characteristics compared in Table 1). Overall, we find limited concordance between the tip rate estimates on our trees and the earlier supertrees of mammals (per-study tip rate correlations of r = 0.60–0.62; Fig 6B). The 221 species identified in the top 1% of tip DR values in our study are similarly recovered in the top percentile at a frequency of 12%, 21%, and 17% for the MRP supertree [87], DNA supertree [92], and consensus timetree [93], respectively. Reducing that comparison to the genus level, those 221 species belong to 46 different genera, of which 22%, 28%, and 41% are similarly recovered with at least one species by those studies (same order). Thus, tip DR estimates from our study best match the DNA supertree at the species level and the consensus timetree at the genus level, although neither are close matches. Relative to each other, the MRP supertree and consensus timetree have the greatest similarity (r = 0.78 in pairwise rates versus r = 0.57–0.59 relative to the DNA supertree). Differences in tree estimation methodologies appear to drive differences in tip rates, although we acknowledge that differences in data availability at the time each study was conducted complicate our comparison.

Overall, the safest statement that we can make is that the tip rate estimates in our mammal trees are substantially different from those of previous trees. However, does that observation equate to our tip rate estimates of mammals being better? Value judgements are difficult in historical biology, in which we lack knowledge of the true evolutionary process. In the absence of simulation studies regarding the efficacy of the backbone-and-patch approach to tree building for recovering true rate dynamics (which would be a welcomed future contribution), we must rely on circumstantial arguments. We can couple the observation of tip DR differences between our trees and previous trees with the following pieces of evidence to make a determination: (1) temporal artifacts are incorporated as a result of supertree merging and polytomies (e.g., rodents Fig 5B); (2) there is a lack of rate uncertainty incorporated in supertrees, particularly when a credible set of trees is not generated (Figs 5 and 6); and (3) the present study was conducted with novel rigor regarding the gathering and cleaning of public DNA sequences, taxonomic reconciliation of synonymous names, supermatrix construction, and use of Bayesian inference methods at levels of the mammalian backbone and subclades (Fig 2).

Departing from past studies while improving the quality of data and inferences argues in favor of our Mammalia phylogenies being able to foster deeper insights into phylogeny-based questions in ecology, evolution, and conservation. Improved understanding of rate dynamics should enable causal hypotheses of biological diversification to be tested with greater reliability [155–157]. Tests that can exclude alternative hypotheses in the presence of realistic tree uncertainty should be viewed as providing durable knowledge regarding the historical trajectory of mammalian evolution.

Limitations

This study was motivated by the clear need for mammalian phylogenetic hypotheses that contain comparably time-scaled branch lengths from root to tip. Until the computational challenges of inferring phylogeny from a genomic matrix of >6,000 species in a joint fossil-calibrated analysis can be overcome and more complete taxon sampling can be obtained, we suggest that the following sources of bias will limit our confidence in the resulting inferences.

Recommended uses of the backbone-and-patch Mammalia trees

We recommend that researchers use the “completed” or “DNA-only” tree sets for addressing questions in which diversification rates or trait evolution are paramount, respectively; when that distinction overlaps (e.g., trait-dependent diversification), we recommend comparing analyses run on tree samples from both sets. In general, all types of analyses should be run on a sample of trees to meaningfully capture uncertainty. Even for questions of ancestral states or character evolution, it is still best to perform topology-based analyses on a sample of DNA-only trees rather than the consensus tree alone [33,34]. Rabosky [175] highlighted that in order to avoid biasing models of character evolution, the unsampled (DNA-missing) species should be avoided because they will move around at random within genus- or family-level constraints [175,176]. Note, however, that species are generally sampled nonrandomly for DNA, so there is an alternative danger of excluding their trait values from analyses. Approaches that apply Rubin’s rules to address missing data in traits and phylogenetic sampling are particularly promising, suggesting that sampling 50–100 trees is sufficient to meaningfully capture parameter uncertainty [177].

Conclusion

The decade plus that has elapsed since the landmark publication of Bininda-Emonds and colleagues [81] highlights the clear need for improved approaches to species-level mammal phylogeny. Our novel, time-calibrated phylogeny incorporating all extant and described species of mammals now enables renewed focus on the causal factors underlying the historical tempo of evolutionary processes. However, the value of continued DNA sequencing for mammal species, as well as the further discovery and cladistic analysis of fossils, should not be understated. Continued improvements to the tree of life, and what we learn about mammalian biodiversity as a result, are directly dependent on the quality of input data. Inferring phylogenies that capture uncertainty in the reconstructed evolutionary process are essential to understanding our mammalian origins.

DNA-gathering pipeline

We used the Basic Local Alignment Search Tool (BLAST) algorithm [105] to query a local copy of NCBI's nucleotide database (downloaded on 20 April 2015), which allowed us to verify standards of homology and orthology among gathered sequence data. Use of BLAST to search for homologous genes avoided name-based searching by taxon or gene and the synonymy issues that entails [178]. We targeted 31 gene fragments commonly sampled among mammals (Table 2), using the family-level supermatrix of Meredith and colleagues [60] as our starting point (22 exons and five noncoding regions). Note that we treated RAG1 as two nonoverlapping regions (RAG1a and RAG1b) to match how researchers have most commonly published these sequences (e.g., GenBank accessions DQ865890 and AY011864 of Didelphis virginiana). To maximize species-level sampling, we targeted four protein-coding mitochondrial genes (mtDNA) in addition to the nuclear genes.

For each gene, we used a set of prevetted sequences, or “baits,” as per-gene BLAST queries of the nucleotide database subset to the NCBI GI list for Mammalia (38,442,994 entries). Given that the per-gene baits were the basis for all downstream work, each bait set was carefully curated to ensure the valid homology of sequences to the targeted gene (e.g., robust alignment, absence of stop codons) while sampling up to one sequence per mammalian family. nDNA baits were taken from the mammal portion of Meredith and colleagues’ [60] DNA alignment, as subset to one family representative per gene. We constructed our own mtDNA baits on a per-gene basis by identifying the longest vetted sequence per family from an output of general NCBI search terms [e.g., for ND2: “((((((txid40674[ORGN] ND2) NOT genome) NOT tRNA) NOT COI) NOT COIII) NOT cytb) NOT d-loop”]. For each gene, we used the “blastn” executable (BLAST+ v.2.2.31) with a coarse E-value of 10 when querying with DNA baits to ensure broad taxonomic coverage. The XML2 output format allowed us to assign NCBI taxonomic information [95] to each resulting hit for subsequent parsing. Using custom Bash scripts, we kept only the unique longest sequence per NCBI taxon ID that was greater than 200 bp in length. Parsing returned a sampling of ≥1 targeted genes for 6,247 unique taxon IDs at the ranks of species and subspecies out of a possible 7,319 such names (85%, as based on the NCBI taxonomy of 20 April 2015). Unaligned FASTA-format files for each gene were then subjected to a taxonomic matchup prior to alignment and further vetting (see S2 Table for steps of successive redundancy reduction, taxonomic matching, and error checking).

This initial procedure yielded direct matches for 4,725 of the 6,247 NCBI names (75%), of which we matched 765 via manual reference to the literature by consulting paper appendices in which the given sequences were published. These manual matches were species (304) or subspecies (410) or had ambiguous epithets (51; “cf.,” “sp.,” etc.). For 135 sequences (77 species), we manually matched accession numbers for which the corresponding NCBI taxon ID pertained to multiple valid species (denoted as “added manually” in S2 Table and S1 Data). Of the 5,490 NCBI names matched directly or manually, there were 1,273 junior synonyms, resulting in a starting list of 4,217 accepted species with ≥1 targeted gene sampled (3,954 species matching the International Union for the Conservation of Nature [IUCN] taxonomy + 263 added species; see below and S2 Fig). These DNA data were the basis for subsequent error checking.

Taxonomy reconciliation and updating

Coherency among species names and their associated data is central to the integrity of any species-level comparative analysis. The NCBI taxonomy associated with our genetic data contained synonymous names and so needed to be vetted against an authoritative list of accepted mammalian species. We chose to initially base this matchup on the IUCN [94] because it (1) followed closely the authority of MSW3 [85], (2) was updated in several cases from MSW3, and (3) was tied to geospatial [179] and species trait [180] resources for downstream analysis. The IUCN base taxonomy contained 5,513 mammal species as downloaded 15 April 2015. We next used a synonym list compiled from several sources (Catalogue of Life [181], MSW3, IUCN; total of 195,562 unique equivalencies; updated from Meyer and colleagues [164]) to match the NCBI species and subspecies names to IUCN.

Our taxonomic matchup also revealed considerable changes to the number of valid mammal species because many were described after the approximately 2004 cutoff date of MSW3 [85] and approximately 2008 cutoff for most of the IUCN list. We made changes to the IUCN base taxonomy as follows: additions of (1) 367 new species, (2) 13 domestic species, and (3) 30 species recently extinct (last approximately 500 years) and subtraction of 12 species synonymized within existing IUCN names. The net change of 398 species resulted in a master taxonomy of 5,911 mammalian species for this study, of which 5,804 species are considered extant (S1 Table and S2 Fig). The Mammal Diversity Database [170,173] (mammaldiversity.org) was an outgrowth of our project, yet it supersedes this total number of species because it continues to update mammalian taxonomy as new literature is published.

Estimation of DNA-only and completed patch clades

Patch clades were estimated with (1) DNA sampling only and no topology constraints (4,098 species, “TopoFree”) or (2) DNA sampling plus taxonomic constraints from the global ML tree to add the remaining unsampled species (5,911 species total, “TopoCons”). For both sets of patch clades, all phylogenies were estimated using identical models of evolution. Parameters of the GTR + G model were estimated independently among the nine partitions, as described above for the global ML tree. Although simpler partitioning strategies might have been selected for some (especially smaller) patch clades, we opted for consistency across all reconstructions. Molecular rate multipliers were estimated for each partition (ratepr = variable) to account for heterotachy (e.g., mtDNA second versus third positions [182,183]). We specified a relaxed clock model of each branch having an independent rate drawn from a gamma distribution [184] and length drawn from a birth–death process (brlenspr = clock:birthdeath). Ultrametric branch lengths were initially estimated in units of expected substitutions per site and subsequently rescaled to absolute time in millions of years from backbone divergence times (see below). Exponential distributions were given to priors for clock rate variance and net diversification (with mean of 0.1) and a beta prior (0–1) on the relative extinction rate. Note that previous backbone-and-patch studies fixed the relative extinction rate to zero, extinctionpr = fixed(0), which resulted in the patch clades being estimated under a pure-birth rather than birth–death process [27,39,40].

We set the in-group sampling probability to the proportion of sampled species per patch for the DNA-only analyses (Table 4) and to 1.0 for the taxonomically completed analyses. For each patch clade, we performed four parallel runs of MrBayes with BEAGLE, each run consisting of four chains of Markov chain Monte Carlo (MCMC; three heated and one cold) and sampled every 10,000 steps for 33,330,000 generations (S1 Text).

For the completed trees, taxonomic constraints for MrBayes were formed with the R package PASTIS [165]. This package reduced the potential for human error while also accounting for nonmonophyletic genera in the global ML tree. Ready-to-execute MrBayes files for each patch clade were generated from the following inputs (see Dryad data: https://doi.org/10.5061/dryad.tb03d03): (1) sequences file of aligned DNA in FASTA format, (2) taxa file of genus membership for all sampled and missing species, (3) missing clades file (if needed) designating where to constrain missing genera, (4) guide tree file giving the relationships of DNA-sampled species (global ML tree, pruned to patch clade species), and (5) template file specifying other MrBayes settings (e.g., rate priors, data partitions). Pruned portions of the global ML tree were not further altered, so all nodes present in that global tree—even those with low BS—were enforced as topology constraints for purposes of adding DNA-missing species to posterior phylogenies. This TopoCons procedure ensured that missing species were added in a phylogenetically informed way, although this has the trade-off of fixing the position of some poorly resolved nodes in the global ML tree. Topologies for those same nodes are sampled probabilistically in the TopoFree DNA-only patch clades.

From the input files, PASTIS adds DNA-missing species to the matrix block with “?” as the character datum for all aligned sites. If left unconstrained, those completed species would be placed at random throughout the posterior sample of trees. The flexible system of “hard,” “negative,” and “partial” constraints in MrBayes, when arranged hierarchically, was used to restrict taxon additions to inside a given clade and outside other clades [50,165]. Partial constraints from PASTIS placed missing species randomly within the least inclusive clade containing all DNA-sampled members of a given genus. For example, the task of constraining 39 missing species of Rattus (of 66 total) was complicated by its paraphyly relative to Diplothrix, Limnomys, Tarsomys, Bandicota, and Nesokia (S3 Data). By using partial constraints, we could constrain missing Rattus species to their paraphyletic grouping across the posterior distribution. Completed species’ branch lengths were drawn from the same birth–death distribution as the rest of the patch clade, biasing PASTIS completions toward rate-constant processes while preserving the taxonomically expected tree shape [27,165]. Therefore, our phylogenies should tend to favor null explanations involving constant-rate species diversification.

Fossil-dated backbone trees

Divergence times and evolutionary relationships among basal lineages of mammals served as the “backbone phylogeny” to which all patch clades were rescaled to absolute time and then joined to form full species-level trees of Mammalia (Fig 2). Fossil information provided the temporal framework for calibrating backbone divergences. Two types of backbones were constructed: (1) ND, using 17 fossil calibrations and one root constraint from Benton and colleagues [139], as augmented by Philips ([185]; S1 Text, section 6 for list of fossil calibrations), and (2) tip-dated (FBD [58]), using the morphological data set of Zhou and colleagues [119] trimmed to 76 fossil (mostly Mesozoic fossils, 66–252 Ma) and 22 extant taxa. In both types of analyses, we focused on a common set of extant taxa to subset the full supermatrix for molecular characters (59 mammals, representing each of the 28 patch clades plus select additional family-level taxa, and one out-group, Anolis carolinensis). Taxa were selected based on their extent of genetic sampling, which for most taxa was >25 of 31 genes (median: 29, range: 3–31). Some additional taxa were selected so that nodes were present for subsequent age constraint in the node-dating analyses or due to their inclusion in the morphological data set of Zhou and colleagues [119].

Node age priors were based on the “best practices for fossil calibrations” recommendations of Parham and colleagues [57], which among other criteria states that (1) fossils should be confidently placed in the crown group of the calibrated node using a formal cladistic analysis of extant and fossil morphological characters and that (2) fossils placed along the stem of a given crown group can only inform the minimum age of the next node back (sister of crown groups [186]). Calibrations were set either as (1) exponential priors (“NDexp”), offset to minimum ages with soft maxima [187] so that the upper 95% of the distribution equaled maximum ages [formula based on the exponential distribution: exp mean = (−ln[0.05]/[max-min])^−1], or (2) as uniform priors (“NDuni”) spanning minima to maxima. These strategies of exponential versus uniform priors were compared to test the sensitivity of dating results. Similar node ages led us to focus on the more conservative NDexp backbone for comparison to the tip-dated backbone and downstream diversification-rate analyses. Analyses were run in MrBayes v.3.2.6 similarly as the patch clades but with the following exceptions: (1) the taxon sampling probability was set to 0.0102 (59 of 5,804 extant species), (2) the sampling strategy was set to “diversity” given our maximization of taxonomic diversity in the backbone phylogeny, and (3) the birth–death clock rate prior was set to a lognormal that assumed each nucleotide site changed one time in the approximately 318 million–year root-to-tip distance (mean = log[1/318], standard deviation = exp[1/318]) [188]. The node age prior was set to “calibrated,” indicating that the probability distribution on terminal and interior node ages was derived from the calibration settings. Note that the birth–death process with variable rates per molecular partition was also implemented here, the same as with the patch clades. We conducted four independent runs of four chains each (three heated and one cold), run for 50,000,000 generations and sampled every 10,000 generations.

For tip dating, we aimed to replicate the study of Zhou and colleagues [119] with the additional “total evidence” perspective of our molecular data set. Tip-dating methods reconstruct the phylogeny of living and fossil taxa according to morphological characters coded for both taxon classes, using the stratigraphic ranges of fossils to inform the birth–death clock model [123,189]. Fossilization can now be parameterized along with birth and death (FBD; [58]), and diversified sampling of extant and fossil taxa can be accommodated in the FBD process [190]. A key contrast in tip dating versus node dating is that fossils impact node ages via their cladistic placement. Thus, different (in some cases fewer) assumptions are involved in tip than in node dating. However, because the former methods are in greater flux [64,121,191], comparing both strategies is useful. Analyses were run in MrBayes v.3.2.6 in a manner analogous to node dating but with the following exceptions: (1) each fossil tip was given a uniform calibration prior between minimum and maximum stratigraphic ages; (2) the “clock:fossilization” branch-length prior was specified, as appropriate for clock trees including fossils [58,190]; and (3) no node calibrations were enforced. Whereas node-dating analyses required topology constraints for each calibration point, none were required in FBD. However, to be sure that our backbone topology matched that recovered in Zhou and colleagues [119], we used hard constraints on the following nodes: in-group from Anolis out-group; crowns (extant taxa) and total groups (with stem fossils) for Placentalia, Marsupialia, and Monotremata; and constraints on the crowns of Theria and Mammalia to ensure the placement of haramiyidans in Mammaliaformes but outside crown mammals ([118,119]; contra to [192]). Note that by exactly following the topology of Zhou and colleagues [119], we constrained two shuotheriids (Shuotherium and Pseudotribos) to be paraphyletic because of an apparent error in their matrix [120]. As with the node-dating runs, for tip dating we conducted four independent runs of four chains each (three heated and one cold), each run for 50,000,000 generations and sampled every 10,000 trees (executables in S5 Data). We estimate that the final runs of node- and tip-dating backbone analyses (not counting troubleshooting) together took 8 weeks in MrBayes using 16 BEAGLE-enabled nodes for a total of approximately 21,500 cpu hours.

Construction of full dated mammalian phylogenies

We summarized the Bayesian patch- and backbone-level runs after discarding the first 25% and 50% of samples as burn-in, respectively. Analyses demonstrated convergent traces in Tracer v.1.5 potential scale reduction factors (PSRFs) of approximately 1 among chains and ESS scores >200 for most parameters, allowing us to combine the four independent runs. Each run was reduced to 2,500 sampled trees after burn-in (from 3,333 for patch clades and 5,000 for backbones), yielding posterior distributions of exactly 10,000 trees upon combining in the Bash shell. For the backbone analyses, fossil-calibrated branch lengths were finalized to absolute time units (millions of years) using the “burntrees.pl” script “--myr” flag (https://github.com/nylander/Burntrees).

To unite the 28 patch clades and the backbone phylogeny (Fig 2A), we first needed to rescale the patch clade distributions of 10,000 trees to absolute time in millions of years. Because patches were estimated in ultrametric but relative units—brlenspr = clock:birthdeath, clockratepr = Fixed(1.0)—this procedure was accomplished with a simple multiplication in R with the “ape” package [193]. We used the rescale-and-graft procedure of previous studies [27,39,188], outlined here:

1. Load all 10,000 trees for a selected backbone analysis and each of the patch clades.
2. Prune all backbones to A. carolinensis and 28 species representatives, one for each of the patch clades.
3. Prune the placeholder taxa from Dermoptera and Platacanthomyidae patch clades.
4. For each of the 10,000 samples of backbones and patches (sequentially),
   1. take one pruned backbone tree; get branching times and identify pairwise MRCA nodes that correspond to the out-group-to-in-group relationship, recording the node age per patch (this root time is used to rescale the patch clade to Ma);
   2. take 28 trees, one from each patch clade; for each, get the relative branching times and divide the age of the in-group crown by that of the root, obtaining the relative scale of the stem edge from out-group to in-group; prune outgroups;
   3. for each patch clade, multiply root time (in Ma) from the backbone by the relative scale of the patch clade to obtain the absolute scale of the stem edge leading to the patch crown; divide all edge lengths by the maximum node height to obtain relative edge lengths of the patch clade and then multiply those by the absolute scale of the stem to rescale all patch clade branches to time in Ma;
   4. use which.edge() in ape to identify the tip edge in the pruned backbone upon which each patch clade is to be grafted; shorten those branches by subtracting the absolute scale of the stem edge per clade; and
   5. use bind.tree() in ape to graft each rescaled patch clade to the corresponding shortened edge of the pruned backbone, thereby forming one uniformly time-scaled tree of 5,911 (or 4,098) species of mammals.
5. Repeat that procedure to construct credible sets of 10,000 phylogenies for each of the backbone analyses (ND and FBD). Note that some authors call these “pseudoposterior” trees because the final trees graft together independent inferences.

Exceptions to the above procedure involved Monotremata, Marsupialia, and Lagomorpha, in which patch clades were rescaled to crown rather than stem divergences because of the use of a distant out-group (Rattus), and some of the bat and rodent patch clades were conditionally assembled because of basal phylogenetic uncertainty (see S1 Text, section 7).

For purposes of visually summarizing the variation in the credible sets of full trees, we (1) constructed MCC consensus trees in TreeAnnotator v.1.8.2 [194] and (2) generated movies to cycle through the variation in tree topology and node ages for trees in each credible set. Mean node ages and 95% HPD intervals were summarized for the backbone-level MCC trees for use in node age comparisons with other studies. For the full species-level trees, we kept the target node heights of the MCC tree rather than annotating node averages (in some cases, node averages would collectively result in negative branch lengths). Although the species-level DNA-only MCC trees are appropriate for some analyses, analyzing a credible set of trees is generally still preferred. The completed MCC trees were for display purposes only (e.g., Fig 1) because the DNA-missing species have variable placements within taxonomic constraints, and thus, no single tree can meaningfully summarize their uncertainty (see Discussion section on “Recommended uses”). Movies of the credible trees’ sets were generated as animated gifs from png-formatted plots of 100 sampled trees using the R package “magick” [195].

Tip-level speciation rates

To characterize species-level variation in speciation rates across Mammalia and compare it among studies, we calculated per-species estimates of expected pure-birth diversification rates for the instantaneous present moment (tips of the tree) using the inverse of the equal splits measure [27,196]. This metric has been called the “DR statistic” and “tip-level diversification rate” because it measures recent diversification processes among extant species [154]. However, to avoid confusion with “net diversification,” for which tip DR is misleading when extinction is very high (relative extinction >0.8 [197]), we here refer to tip DR as a tip-level speciation rate estimator. Tip DR emphasizes geologically recent speciation over deeper-time dynamics, and so it is comparatively less prone to bias from undetected extinction events or nonidentifiability [198] than methods for detecting branch-specific or tree-wide rate shifts [199–201]. We calculate tip DR on full Mammalia phylogenies from the root to each tip as where N_i is the number of edges on the path from tip i to the root, and l_j is the length of edge j. This equation assumes a fully bifurcating tree [24,27]. Because j = 1 is the pendant edge leading to tip i, that branch length carries the greatest weight on the resulting value, with every ensuing rootward edge discounted exponentially as it is shared with other species. Sister species thus have identical tip DR values. Species with the highest tip DR have many short branches shared with other species near the present, implying that recent branching is abundant, whereas low–tip DR species are subtended by long unshared branches (i.e., they are evolutionarily distinct [24]).

We compared our estimates of tip DR calculated across 10,000 trees to those derived from the MRP supertree (1,000 trees; [87]), DNA supertree (1,000 trees; [92]), and consensus timetree of mammals (1 tree; [93]). We summarized the harmonic mean and 95% confidence interval (CI) for the tip DR value of each species in each data set for subsequent plotting. Note that tip DR should only be calculated on “completed” trees that account for all extant and recently extinct species, so we excluded 245 species of Pleistocene extinct mammals in the Faurby and Svenning [92] supertree, retaining only those 97 species presumed to have gone extinct in the last approximately 500 years. Direct matching of species binomial names from our study (n = 5,911) to the other data sets yielded 4,670, 5,329, and 5,033 pairwise comparisons of tip DR, respectively, for these three studies listed above.

Lineage accumulation comparisons

To visualize variation in branching times across Mammalia and among studies, we randomly selected 100 trees from each study and plotted the accumulation of lineages through time from root to tip. We used the ltt.plot() and ltt.lines() functions in the R package ape to output png files for publication.

To assess the congruence of our molecular phylogeny–based rate estimates with the fossil record, we analyzed Mammalia fossil occurrence data from the Paleobiology Database [132], as downloaded on 16 August 2018. Grouping by genus after excluding ichnotaxa and uncertain genera, we recovered 71,928 occurrences of 5,300 genera from the earliest basal Mammaliaformes (e.g., Gondwanadon, late Triassic approximately 235 Ma) to modern genera with fossil records (e.g., Pteropus). The maximum stratigraphic age of the oldest genus per extant mammalian order was used to represent the “fossil stem maximum” age for comparison with our estimates and those of previous studies.

Data availability

All curated data and code are available in the supplementary materials deposited in the Dryad Digital Repository: https://doi.org/10.5061/dryad.tb03d03 [202]. Code for reproducing analyses and figures is also on Github: https://github.com/n8upham/MamPhy_v1. Credible sets of 10,000 trees are available for taxonomic subsetting using the online tool at https://vertlife.org/phylosubsets, and further mammal tree visualizations are presented at http://vertlife.org/data/mammals/.