Morphological data set

Morphological data were obtained from Morphobank project 773 [31]. The unaltered data set contained a total of 4,541 unordered morphological characters that were scored across 86 mammalian taxa, 40 of which are extinct. Specifically, we employed the data set that was used for maximum likelihood analyses in RAxML by O’Leary et al. [31]. This matrix differs from the original matrix by containing no letters for any coded morphological character. This means that the total number of possible character states for any given character was 10 (0–9) instead of 19 (0–9 and: E, F, I, J, L, O, P, Q, X). O’Leary et al. [29] recoded character states that were above 9 as missing (?) while any character states coded as 0–9 were left intact. This change affected twelve morphological characters. This matrix also differs from the original matrix in that all polymorphic character scorings were recoded as missing (?). These changes increased the number of constant characters from 403 to 407.

Each extant order of placental mammals was represented by at least one taxon. Nine of these orders were represented by two or more species. For several of these analyses, the extinct taxa were removed from the matrix, leaving only the 46 extant taxa. Forty-two of these are placentals, two are marsupials, and two are monotremes. Following Beck and Baillie [32] we updated Thomashuxleya with morphological characters from a more recent study [35] and removed 407 characters that were constant across all 86 taxa.

Molecular scaffold

We employed a well-supported molecular scaffold [36] based on Meredith et al.’s [33] phylogenetic analyses of 33,030 base pairs (bp) of protein coding DNA (21 gene segments) and 2,573 bp of DNA from untranslated regions (five gene segments). Clades that were included in our molecular scaffold are also well supported by other molecular studies [16, 37–40]. Our molecular scaffold is illustrated in Fig 1 and included all 46 extant taxa in O’Leary et al.’s [31] data set. This scaffold is fully bifurcating with the exception of the following polytomies: 1. ((Dermoptera + Primates), Scandentia, Glires); 2. (Talpa, Erinaceus, Sorex); 3. ((Saccopteryx + Nycteris), Myotis, Pteronotus); 4. ((Pholidota + Carnivora), Perissodactyla, Cetartiodactyla, Chiroptera); 5. (Hyracoidea, Proboscidea, Sirenia); and 6. (Xenarthra, Afrotheria, Boreoeutheria). We retained these polytomies in all of our scaffold analyses because these nodes were resolved differently in DNA versus amino acid analyses or were not strongly supported (bootstrap support < 90%; Bayesian posterior probabilities < 0.95) [33]. The morphological data of the extant taxa from O’Leary et al. [31] was analyzed with parsimony (see below) and our molecular scaffold to obtain a fully bifurcating molecular scaffold species tree (MSST) (see Fig 1) that was used for ancestral state reconstructions and tree-to-tree comparisons with Robinson and Foulds [41] distances (see below). The molecular scaffold was also employed as a backbone constraint (i.e., containing only a subset of the included taxa) in analyses that included extinct taxa. Finally, extant taxa for each extant placental order were dropped from the molecular scaffold in analyses that treated these orders as pseudoextinct.

Download:

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

Fig 1. Molecular scaffold species tree.

Our fully bifurcating molecular scaffold species tree (20588 steps) that was used for ancestral state reconstructions. The tree was rooted with the monotreme Ornithorhynchus anatinus (platypus). Each internal node has a designated number in red. This number corresponds with the number assigned to each hypothetical predicted ancestor (HPA) in our matrices. These matrices are available in figshare (10.6084/m9.figshare.15183258). Internal nodes that were constrained in the molecular scaffold are denoted with blue circles.

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

Pseudoextinct matrices

O’Leary et al.’s [31] phenomic matrix includes both osteological (1–3,660) and non-osteological (3,661–4,541; total of 881) characters. The latter category includes both soft tissue and behavioral characters. For pseudoextinct taxa, only osteological characters were retained for parsimony analyses. Non-osteological characters were scored as missing (?) for pseudoextinct taxa to mimic the usual absence of these characters in extinct taxa that are only known from fossilized remains. Matrix alterations were performed in Mesquite [42]. After the removal of the 407 constant characters, 4,134 characters remained (3,371 of which were osteological).

Ancestral state reconstructions

The MSST constructed with O’Leary et al.’s [31] morphological data set and our molecular scaffold was loaded into Mesquite. Ancestral state reconstructions (ASR) for each node on our fully bifurcating species tree (Fig 1) were performed with parsimony using the ‘Trace all characters’ command. This procedure was replicated for each pseudoextinct data set, i.e., the order in question was removed from the MSST prior to parsimony ancestral state reconstruction, and ASRs were constructed for a tree that did not include the pseudoextinct order. These reconstructions were exported into Excel where all ambiguously reconstructed characters (i.e., with two or more possible states) were removed and replaced with a ‘?’. Following these alterations, each reconstructed data set was run through the same R script used by Beck and Baillie [32] to ensure that the hypothetical reconstructed ancestors were not coded for any characters that are inapplicable given the reconstructed state at another character. We did not alter this script except to change the name of the input (line 49) and output (line 85) file names. This code is available in figshare (10.6084/m9.figshare.15183372).

We did not explore maximum likelihood ASRs because Beck and Baillie [32] found that the two optimality criterion used in their study (maximum parsimony and maximum likelihood) both fit the Meredith et al. [33] topology equally well when measured by Robinson-Foulds distances [41]. Additionally, in analyses that included all of the hypothetical predicted ancestors, maximum parsimony analyses were more accurate than maximum likelihood analyses in reconstructing monophyletic superordinal clades.

Phylogenetic reconstructions

Each matrix was analyzed in TNT [43] using a traditional search and TBR branch swapping. Morganucodon oehleri was used as the outgroup for all analyses that included fossil taxa. Ornithorhynchus anatinus was used as the outgroup in analyses that omitted fossil taxa. We used maximum parsimony with stepwise addition and 10,000 random taxon addition sequences to search for the shortest tree in each analysis. We did not perform maximum likelihood analyses because there is evidence that parsimony performs better than other optimality criterion for morphological phylogenetic analyses [44–46]. Furthermore, there are concerns about the appropriateness of likelihood models for morphological data since many of these models were designed to fit molecular data [3, 47, 48]. Seven different sets of analyses ((Scaffold + Fossils + HPA); (Scaffold + Fossils); (Scaffold + HPA); (Fossils + HPA); Scaffold; Fossils; HPA) were conducted using different combinations of fossil taxa, HPAs, and our molecular scaffold. In each analysis, all extant taxa were included. Trees were saved as unrooted. For analyses that were performed with the combination of a molecular scaffold + HPAs, we set the maximum number of saved trees in TNT to 100,000 (10,000 random addition sequences x 10 trees per addition sequence). With the exception of Eulipotyphla (99,520 trees), all of the analyses with a molecular scaffold + HPAs resulted in the maximum number of saved trees. The large number of saved trees in the molecular scaffold + HPAs analysis resulted from each HPA having between one and three most parsimonious placements. For example, the HPA for node 65 (Fig 1) between Homo sapiens and Saimiri sciureus has two most parsimonious placements, as a sister taxon to either S. sciureus or (H. sapiens + S. sciureus). However, after pruning all HPAs from these Scaffold + HPA trees, all trees for each pseudoextinct order condensed to a single topology for the 46 extant taxa.

Tree comparisons

Following our analyses, all trees were exported in NEXUS format and imported into PAUP 4.0a [49]. Once in PAUP, any HPAs and/or fossil taxa were pruned from the tree and any duplicate topologies that resulted from this pruning were discarded. All remaining most parsimonious trees, each of which included 46 extant taxa, were then compared to the MSST using Robinson-Foulds (RF) distances [41]. RF distances measure the number of splits (internal branches) that are unique to one of the two trees. All tree comparisons were based on unrooted trees.

The accuracy of morphological phylogenetics

The accurate placement of extinct taxa in the Tree of Life is a fundamental challenge in phylogenetics. Because the placement of most extinct taxa cannot be vetted with molecular data, the accuracy of morphological phylogenetics must be assessed with other approaches such as pseudoextinction analyses that are performed upon extant taxa. Previous pseudoextinction analyses using O’Leary et al.’s [31] morphological data set suggest that morphology alone is insufficient for accurate phylogenetic reconstructions of the majority of pseudoextinct orders [15]. Recently, Beck and Baillie [32] employed a molecular scaffold from Meredith et al. [33] and reconstructed hypothetical predicted ancestors for morphological characters in O’Leary et al.’s [31] data set. These authors concluded that the inclusion of HPAs and fossils in a phylogenetic analysis was sufficient to accurately place all 19 placental orders into their correct superordinal group (Xenarthra, Afrotheria, Laurasiatheria, Euarchontoglires) using morphological data alone. Asher et al. [34] obtained similar results for Glires and found that the inclusion of HPAs nearly tripled the number of well-corroborated groups that were recovered in their morphological phylogenetic analysis of this clade. Importantly, both of these studies estimated and included the immediate HPAs of all extant taxa in their morphological phylogenetic analyses [32, 34]. However, the immediate HPAs of most extinct taxa cannot be reconstructed in scaffold-based analyses because molecular data cannot be procured for these taxa. In this study, we combined HPAs with pseudoextinction analyses to determine if HPAs elsewhere in the tree are sufficient to determine the accurate placement of each ‘extinct’ placental order when its immediate HPA is unknown. We also compared the general performance of molecular scaffolds, fossil taxa, HPAs, and combinations thereof for determining the accurate placement of pseudoextinct taxa.

Our results corroborate the results of Springer et al. [15] and show morphological characters by themselves cannot be relied upon to accurately place placental orders when they are pseudoextinct. Indeed, all seven sets of analyses failed to place the majority of 19 pseudoextinct placental orders next to their correct sister group. Analyses with the molecular scaffold (11 to 12 or 58–63% incorrect sister groups) performed slightly better than analyses that excluded the molecular scaffold (12 to 15 or 63–79% incorrect sister groups). Moreover, at least one of the four superordinal groups was not recovered as monophyletic in 11 to 17 cases (58–89%) when placental orders were treated as pseudoextinct (Table 2). These results stand in sharp contrast to the analysis of Asher et al. [34] that employed HPAs but no pseudoextinction. If pseudoextinct orders are a proxy for extinct orders, the inclusion of HPAs will not result in the accurate placement of pseudoextinct orders when the immediate hypothetical ancestor is unknown. Nevertheless, the inclusion of HPAs, with or without fossils, resulted in greater congruence with well-established clades than analyses with extant taxa only [15] or with the inclusion of fossils but not HPAs (also see [34]).

The challenge of placing insectivores and myrmecophages

Ten placental orders (53%) (Afrosoricida, Carnivora, Chiroptera, Dermoptera, Eulipotyphla, Macroscelidea, Pholidota, Primates, Scandentia, Tubulidentata) were recovered in incorrect phylogenetic positions in all seven analyses. Five of these orders (Afrosoricida, Dermoptera, Eulipotyphla, Macroscelidea, Scandentia) are comprised of taxa that Wagner [50] included in his Insectivora. The polyphyletic origins of Wagner’s Insectivora were recognized by morphologists who removed menotyphlan insectivores from this group and placed them in the orders Dermoptera, Macroscelidea, and Scandentia [51, 52]. However, morphologists continued to support the monophyly of Lipotyphla (solenodons, moles, shrews, hedgehogs, tenrecs, golden moles) [51, 53–55] until this group was shown to be polyphyletic based on molecular data and divided into Afrosoricida and Eulipotyphla [56–59]. Even with the general acceptance of lipotyphlan polyphyly, cladistic analyses of morphology-only data sets have continued to support the monophyly of this clade [31]. Our analyses confirm the difficulty of correctly placing “insectivores” in cladistic analyses with morphological data. Indeed, Afrosoricida and Eulipotyphla always join each other when one of these groups is pseudoextinct.

Beyond the various “insectivore” groups (Insectivora, Menotyphla, Lipotyphla) that have subsequently been shown to be polyphyletic, Edentata is a polyphyletic taxon that is historically comprised of the myrmecophagous orders Xenarthra, Pholidota [22–24] and sometimes Tubulidentata [31]. Two of these orders (Pholidota, Tubulidentata) are among those that were always misplaced in our pseudoextinction analyses. Pholidotans were always the sister taxon to Xenarthra or Pilosa when pseudoextinct. Tubulidentata, in turn, were usually the sister taxon to other myrmecophagous taxa (i.e., Pholidota, Xenarthra, Pholidota + Xenarthra) although in two cases (Scaffold, Scaffold + HPAs) were the sister taxon to Euungulata (Cetartiodactyla + Perissodactyla). With the exception of Tubulidentata + Euungulata, these results represent different incarnations of Edentata and underscore the challenges that are associated with accurately placing myrmecophagous taxa in morphological phylogenetic analyses. The order Carnivora is represented by Canis familiaris (dog) and Felis silvestris (European wildcat) in O’Leary et al.’s [29] data set. With a few exceptions (e.g., Proteles cristata), the members of Carnivora are not insectivorous/myrmecophagous. However, the sister taxon to Carnivora is an edentulous, myrmecophagous order (Pholidota) [33] and was misplaced in all seven analyses when pseudoextinct. This result further highlights the difficulty of reconstructing clades such as Ferae (Carnivora + Pholidota) that include myrmecophagous orders.

By contrast with our results for Pholidota and Tubulidentata, the myrmecophagous Tamandua tetradactyla (Order Pilosa) and largely myrmecophagous Dasypus novemcinctus (Order Cingulata) were often recovered as sister groups to each other when pseudoextinct (Cingulata in 6 of 7 analyses; Pilosa in 5 of 7 analyses; see Table 2). Previously Springer et al. [15] performed a pseudoextinction analysis with O’Leary et al.’s [31] data set and treated Cingulata and Pilosa as a single order (Xenarthra). Pseudoextinction analyses with Xenarthra resulted in this taxon jumping inside of Afrotheria and joining with Tubulidentata [15]. Similarly, Springer et al. [12] performed pseudoextinction analyses with Xenarthra and a much smaller data set (185 osteological characters; [60]), and recovered Cingulata as sister to Tubulidentata and Pilosa as sister to Pholidota. These results suggest that morphological data will recover more accurate positions for extinct taxa when a suitably close living relative is also included in the analysis.

Effect of HPAs

The inclusion of the HPAs resulted in a greater number of pseudoextinct orders that were recovered as monophyletic. For two analyses, Scaffold + Fossils + HPA and Fossils + HPA, all pseudoextinct orders were recovered as monophyletic. In the Scaffold + HPA and HPA only analyses only one pseudoextinct order (5%) (Cetartiodactyla) was not monophyletic. In the HPA only analyses Cetacea was recovered outside of an ungulate clade that included Perissodactyla and terrestrial Cetartiodactyla. In the Scaffold + HPA analyses, Cetacea was recovered as the sister group to Sirenia. The increase in the number of monophyletic (pseudoextinct) orders that resulted from the inclusion of HPAs suggests that increased fossil sampling may also increase the accuracy of morphological phylogenetic analyses if a sufficient number of transitional fossils with appropriate combinations of primitive and derived characters are included in the analysis.

However, HPAs may not replicate extinct taxa in the fossil record [34]. Thus, new fossil discoveries may not have the same positive effect as HPAs when they are included in a cladistic analysis. Additionally, the fossil record is highly incomplete and many gaps in the fossil record will likely be permanent due to environmental factors that did not favor fossilization [61]. Thus, even with additional fossil discoveries, it is unlikely that all ancestors for a large group will ever be known. Our results show that the absence of even a single ancestor can be severely detrimental to the accurate phylogenetic reconstruction of some pseudoextinct groups. We acknowledge that the analyses of Beck and Baillie [32] represent a best-case scenario for future fossil discoveries. It remains to be tested if close relatives that are not ancestral would be as beneficial as ancestral species to improve the performance of morphological phylogenetic analyses.