The Interrelationships of Placental Mammals and the Limits of Phylogenetic Inference

1The quest for the root of placental mammal phylogeny has achieved the status of an iconic controversy, with three principal competing hypotheses that resolve either 1) Xenarthra (e.g., armadillos and sloths;;; ), 2) Afrotheria (e.g., elephants and tenrecs;;;;;; ), or 3) Atlantogenata (i.e., Xenarthra plus Afrotheria;;;;;; ) as the sister to all other placentals (fig. 1). Previous analyses have found support for all three hypotheses, leading some to conclude that this phylogenetic problem is impossible to resolve. This has been considered a consequence of incomplete lineage sorting (ILS;;;; ), reflected in large scale gene tree heterogeneity, a result of the apparent rapidity of successive vicariance-driven divergence events associated with the fragmentation of the Pangaean and Gondwanan supercontinents. Thus, if placental mammals evolved extremely rapidly, then the root of the placental radiation may be theoretically unresolvable, as it was not strictly bifurcating in the first instance. However, it is possible that phylogenetic resolution has been precluded by practical constraints, which include the availability of adequate models of molecular evolution, compositional biases, and/or long branch attraction, and computational limitations on the scale of molecular sequence data sets with limited gene and/or taxon sampling. Resolution among these three competing hypotheses is essential to understand the evolutionary origin and diversification of placentals, the most phenotypically diverse group of vertebrates, occupying terrestrial, aerial, and aquatic ecological niches, with body sizes spanning several orders of magnitude and which were accompanied by both large scale genomic (e.g., transposable elements,; conserved noncoding RNAs, ) and morphological (e.g., the placenta; ) innovation.

2In an attempt to resolve this phylogenetic controversy, we undertook analyses of two genome-scale data sets representing both coding and noncoding regions of the genome: a 21.4 million nucleotide superalignment of 14,631 genes from 36 taxa, and a 16,050 nucleotide superalignment of 239 pre-miRNAs from 39 taxa. In addition, we reanalyzed the data from three recent analyses that obtained results incongruent with those from our protein coding and nonprotein coding data sets, and tested the extent to which morphological data can inform mammal phylogenetics using the 4,541 character data set of.

3We performed phylogenetic analyses of two nucleotide data sets and three amino acid data sets. The nucleotide data sets were a superalignment of 14,631 protein-coding genes and 36 taxa (totaling 32,116,455), and a superalignment of pre-miRNA sequences comprising 16,050 sites and 42 taxa. The three amino acid data sets were the 11,365 amino acid data set of, the AT-rich amino acid data set of, and the amino acid data set of. For all considered data sets Posterior Predictive Analysis (PPA) of biochemical specificity was performed to investigate whether standard, compositionally site-homogeneous, models (e.g., general time reversible [GTR] and Whelan and Goldman [WAG]) provided an adequate fit to the data or whether a more complex (compositionally site–heterogeneous) model (e.g., CAT–GTR;; ) was necessary to adequately fit the data. For the nucleotide and microRNA (miRNA) data sets two models were tested, the GTR+G model and CAT–GTR+G. For the amino acid data sets PPA was used to compare the model used in the original studies (Jones, Taylor, and Thorton [JTT]+G; LG+G; and WAG+G), against the CAT–GTR+G model. PPA was performed using the serial version of Phylobayes 3.3f (following suggestions from Nicolas Lartillot) using data sets that were subsampled to include a set of approximately 5,000 randomly selected characters. The final number of characters is variable (but comparable) across the different data sets, because of the subsampling strategy we used. However, this is not important as models are compared on the same data sets and not across data sets.

4Total data set size and percentage of missing data is record in table 1. All Bayesian analyses were performed using the CAT–GTR+G model and implemented with the MPI version of the software Phylobayes (Phylobayes MPI 1.5a; ). For all Phylobayes analyses two chains were run. Burn-in varied and all chains were run until convergence (which was tested using the BPCOMP software, which is part of the Phylobayes suite). Following the Phylobayes manual, chains were considered to have converged on the same solution when the Maxdiff (maximal difference between observed bipartitions) dropped below 0.2. Maximum Likelihood analyses were performed using RAxML under a GTR+G model, and the bootstrap (100 replicates) was used to estimate support.

5The genome alignment of, comprising 36 taxa and 14,631 protein-coding genes was used. Codon sequences were aligned using PRANK with no guide tree to minimize bias associated with any guide tree, although we note that alternate alignment software will generate alternate alignments and subsequent analyses should examine whether such alignments affect our results. The first and second codon positions of all genes were concatenated into a single partition (21,410,970 nt). Because of computational limitations, the full data set could only be analyzed using maximum likelihood. We investigated whether the results of our maximum likelihood GTR+G analyses were supported also under CAT–GTR+G, but because a CAT–GTR+G analysis of the entire superalignment is unfeasible, we removed all the constant and parsimony uninformative sites prior to the analysis. We recognize that this is not ideal, as it can introduce biases and this analysis can consequently be considered to have only an exploratory nature. Initial CAT–GTR+G analyses included all the taxa but did not converge. Inspection of the two chains showed that the horse and tree shrew were unstable within Boreoeutheria. As these taxa are irrelevant to investigate the relationships at the root of the placental tree (Boreoeutheria was monophyletic in both chain and with a posterior probability of 1), we repeated analyses excluding these two taxa. This analysis converged on the same topology within 150 generations (with a Maximal Difference between observed bipartitions dropping to zero).

6After having run our phylogenetic analyses we investigated whether the data could significantly discriminate between alternative hypotheses of placental relationships. As CAT–GTR+G and GTR+G supported the same tree for the nucleotide data set these analyses were only explicitly performed under maximum likelihood using the GTR+G model. To do so, the three competing hypotheses were fixed and compared using the approximately unbiased (AU)-Test. Sitewise likelihood values were obtained (under each considered hypothesis of placental relationships) using BASEML, and CONSEL was used to calculate the AU test. Because of computational limitations AU tests was only performed using the superalignment, and not on the 14,631 individual gene alignments constituting our superalignment. For the gene-by-gene analyses a reduced data set of 11,169 genes was used so that every gene had at least one non-placental outgroup, a Xenarthran, Atlantogenatan, and Boreoeutherian present in the alignment so that the tree could not only be rooted but was also informative as to the relationships between these key clades. For each gene we then estimated the likelihood of each considered tree and performed two different analyses. First, we calculated how many genes supported each alternative hypothesis without considering whether the differences in likelihood between compared trees have been significant. This identified the number of genes for which each considered topology is optimal. Subsequently the Akaike Information Criterion (AIC) test was used to determine whether the genes supporting each specific tree topology, supported that topolgy significantly better than the other tree topologies.

7The reduced data set of 11,169 genes from the gene-by-gene analyses (see above) was used to define the set of unbinned gene trees. We also used a statistical binning pipeline with support threshold set to 50% to create 2,513 bins of genes (1,373 bins with four genes, 1,139 bins with five genes, and one bin with six genes) and estimated a supergene tree for each bin. ASTRAL version 4.7.6 was run on both sets of inputs: the 11,169 unbinned gene trees, and the 2,513 supergene trees, weighting each supergene tree by size of the corresponding bin (weighted statistical binning;; ). To test for the number of gene trees that supported each hypothesis with support above 50% or 75% threshold, we first restricted each gene tree to branches that have support above the chosen threshold. We then compared each collapsed gene tree against three unresolved trees that represented the three hypotheses. A gene tree can either reject all three hypotheses (i.e., when Xenarthra, Afrotheria, Boreoeutheria, or the branch uniting the three outgroups are rejected), or be indecisive (i.e., be compatible with all three hypotheses; this happens when in the collapsed gene tree, the relationship between Xenarthra, Afrotheria, Boreoeutheria is unresolved), or can support one of the three hypotheses. Thus, five outcomes are possible, and we note the percentage of times each outcome is observed. We also note the percentage of gene trees that support each of the three hypotheses out of those that support just one hypothesis. This produces three estimated probabilities, one for each hypothesis, and we can convert these probabilities to coalescent unit branch lengths by calculating −ln (3/2* (1−p)) where p is the probability of a hypothesis. For example, for unbinned gene trees, out of 3,495 genes that exclusively supported one the three hypotheses with at least 50% BS, 48.4% of them supported Atlantogenata, which puts the branch length in coalescent units at −ln (3/2* (1–0.484)) = 0.257. Using 75% threshold with unbinned gene trees results in a length of 0.415, and using supergene trees with 50% and 75% threshold result in lengths of 0.135 and 0.192, respectively.

8Small RNA libraries were generated from whole juvenile specimens of Armadillo (Dasypus novemcinctus), Rabbit (Oryctolagus cuniculus), and Guniea Pig (Cavia porcellus) using the Illumina Tru-seq small RNA prep kits. In brief, this process involves taking 1 µg of total RNA and adding 5′- and 3′-adapters, which were then reverse transcribed, barcoded, and amplified using polymerase chain reaction. The sample was run out on a Novex 6% TBE Page gel using electrophoresis allowing size fractionation of the sample. The relevant size fraction will be excised and eluted overnight to increase total product. The eluate will be precipitated using EtOH, glycogen, and sodium acetate for 24 h before being resuspended and submitted for sequencing on a GAIIx sequencer at the University if Bristol Transcriptomics Facility. Total read counts were approximately 22M for Armadillo, approximately 13.5M for Guinea Pig, and approximately 21M for Rabbit, and the data processed using in-house algorithms. These read data were used to verify the mature and star reads and hence the end of the pre sequence, which was used for the pre-mir alignments and have been deposited in miRBase. In addition, BLAST searches were conducted for an additional 42 taxa to identify additional miRNA loci. Orthology for each individual miRNA was checked using both distance and, when possible, syntenic analysis. Each individual pre-miRNA from the 42 taxa analyzed was concatenated into the tetrapod superalignment of and analyzed as a standard superalignment comprising 15,590 sites and 42 taxa, using the GTR+G model.

9Several recent studies addressed the relationships among the placental mammals finding contradictory results. A feature characterizing these studies is the heterogeneity in the choice of the model used for phylogenetic analyses, and the fact that in all cases the substitution model used to analyzed the data was selected in either a subjective way or from a subset of models that did not include well-performing (parameter rich) site-heterogeneous models. Following the results of our PPA (see above), which showed that the models used in the original studies did not fit the data adequately, the three data sets associated with these studies (the 11,365 amino acid data set of, the AT-rich amino acid data set of ], and the amino acid data set of ) were reanalyzed under the site-heterogeneous CAT–GTR+G model.

10recently presented a 4,541 character morphological data set. We tested whether this morphological data set could distinguish between the three alternative hypotheses of placental relationships. As in the case of the nucleotide data set the AU-Test was used (implemented in CONSEL), with character-wise likelihood values estimated in RaXML under the MKv model.

11The 21m nucletotide alignment was used for the molecular clock analysis. This alignment has previously been used, however, the discovery of new fossil material, as well as revised stratigraphy and phylogenetic placement of taxa means that 20 of the 23 calibration points shared between studies had to be revised (table 2). The previously unpublished calibration on node 37 is justified below following best practice guidelines.

12Fossil Taxon and Specimen: Haramiyavia clemmenseni (Museum of Comparative Zoology MCZ 7/G95) from the Tait Bjerg Beds, Ørsted Dal Member of the Fleming Fjord Formation with an age corresponding to the Late Triassic (?Norian-Rhaetic).

13Phylogenetic Justification: Prior to the discovery of Haramiyavia clemmenseni, haramiyids were known from two genera. However, the taxonomic status of these genera was uncertain, and while H. clemmenseni exhibited highly specialized dentition it also retained features of the jaw and post-dentary apparatus that indicated a position among stem mammals, cladistically more basal than crown Mammalia, i.e., the clade encompassing monotremes and therians. Some recent phylogenetic studies have placed Haramiyavia as sister taxon to multituberculates, which are closer to therians than to monotremes and thereby within crown Mammalia. In contrast, other studies argue that the anatomical similarities between haramiyids and multituberculates are convergent. We tentatively use Triassic haramiyids as a minimum calibration for Mammalia but are keen to see future, more thorough phylogenetic tests of haramiyid affinities.

14Minimum Age: 201.1 Ma

15Soft Maximum Age: 252.23 Ma

16Age Justification: At present Haramiyavia clemmenseni is the oldest known haramiyid from the Tait Bjerg Beds, Ørsted Dal Member of the Fleming Fjord Formation with an age corresponding to the Late Triassic (?Norian-Rhaetic). This stage (Rhaetic) has a minimum bound of 201.3 Ma ± 0.2 Myr and so the soft minima is set at 201.1 Ma.

17Broader Justification: Hadrocodium and Docodonta are the closest relatives to crown mammals. Hadrocodium is from the early Jurassic of Yunnan Province, China (Sinemurian; ), and the oldest docodonts are from the Bathonian of Europe, both of which are younger than Haramiyavia. More distantly related taxa such as Morganucodontidae, Sinoconodon, and Adelobasileus, are known from the late Triassic and early Jurassic and are contemporaneous with Haramiyavia, implying substantial ghost lineages in many of these taxa, as such a broad prior is used, setting the soft maxima at the PT extinction, dated at the base of the Induan, 252.17 Ma ± 0.06 Myr and so the soft maxima is set at 252.23 Ma.

18The molecular clock analysis was conducted with MCMCTREE v. 4.8 a, using the approximate likelihood method by calculating the maximum-likelihood estimates of the branch lengths, the gradient vector and Hessian matrix, using BASEML, under the HKY+G4 model. We then used the Markov chain Monte Carlo algorithm to estimate divergence times on the constrained tree topology with two separate runs being performed. The auto-correlated rates model was used to specify the prior of rates, and we followed for other parameters, that is; the time unit was 100 Myr; a diffuse gamma prior G(1, 1) was used for the overall substitution rate; a rate drift parameter σ2 was assigned G(1, 1); and the parameters of the birth–death process with species sampling were fixed at λ = μ = 1 and ρ = 0. The alignment was analyzed as a single partition and we conducted 2,000,000 iterations, sampling every 200 a burn-in of 25%, and with both runs being concatenated post burn-in, after thinning down to 10,000 samples per run, to provide the final posterior values.

19A fully resolved phylogeny with 100% support for both a sister group relationship between Afrotheria and Xenarthra (Atlantogenata) and between Atlantogenata and Boreoeutheria (fig. 2, left; supplementary fig. S1, Supplementary Material online) was recovered in the analysis of the 21.4 million nucleotide alignment (first and second nucleotide positions) using a single GTR+G model. Of the 35 internal nodes, 32 were recovered with 100% support. Further analyses were performed using the compositionally site–heterogeneous CAT–GTR+G model, which accommodates among-site amino acid (and nucleotide) compositional heterogeneity. This analysis recovered the same topology with all nodes exhibiting 100% support (fig. 2, left; supplementary fig. S2, Supplementary Material online). Unambiguous statistical support for Atlantogenata was confirmed using the AU test, which assesses the level of support for each topology through a site-by-site analysis of the entire data set. The results of this analysis rejected basal positions for both Afrotheria and Xenarthra (P ≤ 0.01) in favor of Atlantogenata (P ≥ 0.99) (table 3).

20Despite strong support for Atlantogenata, we decided to investigate the level of support for each of the three topologies from the individual genes. We therefore removed all of those genes that were unique to individual lineages, that is, Euarchontoglires, Laurasiatheria, Primates etc., or that were not represented by at least one member of Xenarthra, Atlantogenata, Boreoeutheria, and a nonplacental outgroup. This was done so that each individual gene had the potential as to be informative to the placental root, resulting in a reduced data set of 11,169 genes. The number of individual gene trees recovering alternative topologies (albeit not necessarily with high support) is comparable: Atlantogenata (∼33.88%), Afrotheria (∼33.84%), Xenarthra (∼29.9%), and indecisive (∼2.3%) (see table 3). These results could be interpreted to support the prevailing view that the early phylogenetic history of placental mammals was such a rapid radiation that it was not strictly bifurcating. However, 99.4% of the genes fail to discriminate among the competing hypotheses with statistical significance as measured by the AIC test, leaving only 0.2% of genes supporting Atlantogenata, 0.12% supporting Afrotheria, and 0.22% supporting Xenarthra (table 3). Thus, the distribution of support for competing topologies largely reflects the weak phylogenetic signal present in any single gene alignment, rather than suggesting a hard polytomy or very high levels of ILS.

21It is known that concatenation analyses, such as those performed here, can be statistically inconsistent or even positively misleading in the presence of sufficient levels of ILS. Thus, we further tested the robustness of our phylogeny through the use of ASTRAL-2, a coalescent-based species tree estimation method that is robust to the presence of ILS. We also explored the use of weighted statistical binning, a technique designed to improve species tree estimation when gene trees have poor resolution. Thus, we used ASTRAL with and without weighted statistical binning, applied to the same 11,169 genes used in the gene-by-gene analysis described earlier.

22In both cases a fully resolved tree with 100% support for Atlantogenata (fig. 2, left; supplementary fig. S3, Supplementary Material online) was returned, supporting the concatenation analysis. After restricting analyses to the set of gene trees with high bootstrap support (50% or 75%) for one of the considered hypotheses, support for Atlantogenata was strengthened (supplementary fig. S4, Supplementary Material online). For example, 48% of the unbinned genes and 42% of the binned supergenes that met the 50% bootstrap support threshold supported Atlantogenata, with almost equal numbers of genes supporting an Afrotheran (26% or 30%) or Xenarthran outgroup (26% or 28%). When the level of bootstrap support necessary for the gene trees to be included in the analyses was increased to 75%, the preference for Atlantogenata further increased to 56% of the unbinned genes, and 45% of the binned supergenes, with corresponding decreases in the levels of support for Afrotheria or Xenarthra (fig. 3). This suggests that some (and perhaps much) of the incongruence observed across the gene trees is the result of stochastic errors in gene tree estimation, not ILS. When restricted to gene tree branches that have bootstrap support above 50% or 75%, the branch length for the Atlantogenata group is between 0.14 and 0.42 coalescent units (depending on the threshold and/or the type of gene trees used; see table 4 and Materials and Methods). Critically, the highest levels of support and longest branch lengths in terms of coalescent units for Atlantogenata are returned when we analyze the data using unbinned gene trees. Our estimated coalescent unit branch lengths point to a short branch, but not an extremely short branch that would violate the hypotheses of a strictly bifurcating tree. These results are largely congruent with concatenation analyses, and suggest that the amount of discordance due to ILS is not sufficient to mislead the concatenation analysis. Thus, although the two analyses are based on data sets of different sizes (11k and 14k genes, respectively), both types of analysis—coalescent-based and concatenation—are highly congruent, and both provide high support for Atlantogenata.

23In addition to protein-coding genes, we also assembled a concatenated superalignment of 239 noncoding RNA miRNAs consisting of 16,050 nt, which was analyzed under the GTR+G model (see table 5). This miRNA data set provides a second independent molecular data set, that of noncoding RNA genes, to complement protein-coding gene analyses, and these data can be analyzed using the same model based approaches. Such an approach has been shown previously to be suitable in resolving interspecies relationships among reptiles, primates, nematodes, and drosophilids. Our pre-miRNA superalignment recovered a fully resolved tree with an Atlantogenata outgroup exhibiting a posterior probability of 0.97 (fig. 2 right; supplementary fig. S5, Supplementary Material online), in agreement with the protein-coding gene analyses. We again used an AU test to investigate site-by-site support for the three topologies on the entire data set with the results significantly rejecting Afrotheria with P = 0.028 (Xenarthra P = 0.250; Atlantogenata P = 0.795), once more providing support against a hard polytomy.

24Given the consistent support in our two data sets for Atlantogenata, we explored why some previous data sets did not find support for this rooting. Amino acid data sets have yielded support for Afrotheria and Xenarthra, and analysis of an AT-rich amino acid data set supported Afrotheria. We focused on model selection and using PPA we showed that the models used in the original studies (WAG2000+G, JTT+G, LG+G, respectively) did not adequately fit the data (see table 5). In contrast, for each of these three data sets, the compositionally site-heterogeneous CAT–GTR+G model was found to be a satisfactory fit to the data. Reanalysis of all three data sets using the CAT–GTR+G model found variable support for Atlantogenata (supplementary figs. S6–S8, Supplementary Material online), and not for the relationships reported in the original studies, undermining their conclusions. Support values for an Atlantogenata root vary considerably between the three reanalyses with values of 1 (; supplementary fig. S6, Supplementary Material online), 0.79 (; supplementary fig. S7, Supplementary Material online), and 0.5 (; supplementary fig. S8, Supplementary Material online). While a support value of 50% is uninformative the original paper had a bootstrap support of 100% for Afrotheria. Thus, although this reanalysis does not have high support the use of a better fitting model fundamentally overturned the previous hypothesis, which was itself very highly supported. Likewise, the results of, which previously supported Xenarthra, were overturned to support Atlantogenata. Furthermore, these two data sets with the lowest levels of support either contained low numbers of loci (27 nuclear genes) as in or sampled a nonrandom selection of genes, focusing on AT-rich genes as in, such approaches are likely to exacerbate phylogenetic artefacts through both compositional and long branch attraction.

25In addition to their amino acid data set, also used a 4,541 character morphological datamatrix. When this matrix was analyzed using the AU test in RAxML (with a constraint tree to make Afrotheria, Xenarthra, and Atlantogenata monophyletic) the morphological data set was unable to distinguish between the three competing hypotheses (Afrotheria P = 0.288, Xenarthra P = 0.212, and Atlantogenata P = 0.363). Thus, when analyzed in isolation, the morphological data are indecisive concerning the earliest diverging lineage of placental mammals.

26We estimate the mean divergence times for crown Theria as 164 Ma (CI = 157–170 Ma), crown Placentalia as 93 Ma (CI = 86–100 Ma), and crown Atlantogenata as 90 Ma (CI = 84–97 Ma) (fig. 4 and table 6). These dates are considerably younger than some studies, older than others, and congruent with others still. As expected, our revised calibrations, older than those employed by ), have the effect of making the posterior ages slightly older (Placentalia and Atlantogenata increase in mean age by 3.1 and 2.5 Myr, respectively), while the 95% CI broadens from 88.3–91.6 to 86.5–99.9 Ma in placentals, and 85.9–89.1 to 83.7–96.5 Ma in Atlantogenata. This broadening in the 95% CI reflects the use of a single data partition, in comparison to in which 20 partitions were used. We estimate diversification of placental orders overlapping the K-Pg mass extinction event at 66 Ma, with all placental orders diversifying between 76 and 51 Ma.

27We advocate the use of providing details of not only the combined posterior, but also the marginal prior, which is an analysis run without the sequence data so that the effect of all calibrations can be assessed, as the marginal prior for any node can differ from the original fossil calibration. Here we observe that the marginal prior closely approximates (<2 Myr) the calibration points at the majority of nodes, with only four exceptions. The only substantial deviation (>4 Myr) is with the soft minima on node 56 (Euarchontoglires). Such results show that the priors were performing as expected based upon the initial fossil calibrations.

28Thus, far from an intractable phylogenetic problem, it is evident that conflicting placental phylogenies have been a consequence of the use of poorly fitting evolutionary models. Evidently, there was some gene tree heterogeneity caused by ILS during placental diversification. However, we can reject the view that this was so rampant as to obscure the fundamental relationships among placental mammals. Instead, our results demonstrate that the primary evidence on which such ideas are based, that is, an equal number of genes supporting mutually exclusive topologies, is the consequence of weak signal in single gene alignments rather than the result of ILS alone. As articulated elsewhere (e.g.,;; ), isolated, historical signal becomes stronger when individual partitions (such as gene alignments) are combined. Thus, we reject the view that the root of the placental mammal tree is an unresolvable polytomy, concluding instead that it is correctly resolved as a fundamental divergence between Atlantogenata and Boreoeutheria.

29We do not doubt that evidence of ILS reflects the fact that the initial diversification of placentals was rapid as is observed in our molecular clock analysis, the results of which are comparable to those reported elsewhere. Although the discovery of several recent fossils has led to the calibrations being revised substantially with the minimum ages for the root (Mammalia) and Theria being pushed back 38.2 and 32.3 Myr, respectively, whereas the maxima age for Placentalia was pushed back by 33.1 Myr. Yet, such revisions had only minor changes in the estimated mean age of diversification for Placentalia (+3.1 Myr), Atlantogenata (+2.5 Myr), and Boreoeutheria (+2.6 Myr), however, larger changes were observed for the Mammalia (+38.9 Myr), and Theria (−11.5 Myr) in comparison to the results of. These dates support dispersal, rather than vicariance, as the underlying mechanism in placental mammal biogeography as they postdate not only the fragmentation of Pangaea, but also the later splitting of Gondwana due to the opening of the Atlantic Ocean.

30Previous studies have suggested a clear pattern of biogeographic diversification for placentals into four principle lineages (Afrotheria, Xenarthra, Laurasiatheria, and Euarchontoglires) caused by drift-vicariance, which followed the continental breakup of Pangaea into the northern continent of Laurasia (Laurasiatheria + Euarchontoglires) and a Southern Gondwanan continent (Afrotheria and Xenarthra) in the Jurassic (201.3–145 Ma). This was followed by the later breakup of Gondwana into South America (Xenarthra) and Africa (Afrotheria) due to the opening of the Atlantic during the Cretaceous approximately 110 Ma. Recent analyses of global plate tectonics suggests these dates for the complete breakup of Gondwana into S. America and Africa are too old and that this separation was fully complete by 100; ). However, these dates not only predate our mean divergence time for the divergence of Afrotheria from Xenarthra by approximately 10 Myr, but they also lie outside of the 95% HPD (83.73–96.54 Ma), suggesting dispersal by a group of stem Xenarthrans across the Atlantic. While dispersal across the proto Atlantic Ocean may seem unpalatable, the scale of the Atlantic ocean barrier in the Late Cretaceous (fig. 4) was far less significant than that between Africa and Madagascar which has, nevertheless, witnessed multiple post-Mesozoic dispersal events of placentals, including tenrecs, rodents, primates, and carnivores. Oceanic dispersal of rodents and primates across the South Atlantic during the Eocene (when the overwater distance between Africa and S. America was wider compared with the Cretaceous) is also uncontroversial.

31With the resolution of the evolutionary relationships among Afrotheria, Boreoutheria and Xenarthra, attention must now turn to resolving the problematic relationships within Laurasiatheria and to understanding of the role of dispersal in effecting placental diversification. The results of both our RAxML and ASTRAL analyses as well as the reanalyses of and place the tree shrew as sister taxa to Glires, and the horse in an Euungulata clade, and suggests that classical groupings such as Euarchonta and Ferungulata are not supported. Although such results have been presented before this is an area of significant conflict between previously published studies. It is these two rogue taxa (tree shrew and horse) which are the cause of alternate tree topologies, and it is no surprise that the same two taxa were the ones that needed to be removed from our phylobayes analysis as they prevented the runs from converging. In future increased taxonomic sampling of additional perissodactyl lineages, that is, Equidae (donkeys, and zebras), Rhinocerotidae (rhinos), and Tapiridae (tapirs) as well as Scandentia lineages, that is, Anathana (Madras treeshrew), Dendrogale (Bornean smooth-tailed treeshrew), and Ptilocercus (Pen-tailed treeshrew), will lead to increased confidence in the phylogenetic placement of these lineages. While a better understanding for the role of dispersal through not only the late Mesozoic but also the Paleogene (or early Cenozoic) can be achieved through a more precise understanding of the geography including sea-level changes, and not merely the tectonics and biogeography through this interval. In addition the inclusion of fossils within analyses of their living relatives needs to become more widespread, allowing not only greater precision in divergence time estimation through the use of tip dating in molecular clock analyses, but also to better understand the pattern of character acquisition, and changes in diversity, either to identify diversification rate shifts or broader patterns of biological diversity.

32The results of our study suggest that other seemingly intractable phylogenetic debates, such as the position of ctenophores, chaetognaths, Acoelomorpha, and the relationships among lophotrochozoans, may be solvable by combining genome-scale data sets with realistic models of molecular evolution and rigorous coalescent-based species tree estimation methods.

33Supplementary figures S1–S8 are available at Genome Biology and Evolution online (http://www.gbe.oxfordjournals.org/).