Introduction
Phylogenetic methods that employ parametric evolutionary models, such as the Lewis Mk model (Lewis Reference Lewis2001), have recently begun to emerge as important tools for addressing paleobiological issues at macroevolutionary timescales (Wright and Hillis Reference Wright and Hillis2014; Puttick et al. Reference Puttick, O'Reilly, Tanner, Fleming, Clark, Holloway, Lozano-Fernandez, Parry, Tarver, Pisani and Donoghue2017). Probabilistic approaches appear promising in their performance relative to earlier cladistic methods and offer several important advantages. These include their tendency to increase the clarity with which models are specified and their ability to weigh evidence under competing evolutionary models using modern inferential machinery. While probabilistic approaches have long been dominant in molecular phylogenetics, their application to morphological data sets remains immature. One topic that has only recently been addressed in parametric approaches to phylogenetic inference has been the incorporation of paleobiological evolutionary patterns, such as budding cladogenesis and anagenesis.
Traditional approaches usually depict evolutionary relationships in terms of branching, or in paleobiological terms, splitting cladogenesis. However, when analyzing the fossil record, it is often desirable to entertain hypotheses wherein taxa occurring earlier in time are directly ancestral to taxa that occur later in time. The presence of such directly ancestral taxa is expected in the fossil record (Foote Reference Foote1996) and is central to identification of paleobiological patterns such as budding cladogenesis and anagenesis. Characterization of evolutionary mode by distinguishing splitting cladogenesis from budding and anagenesis in the fossil record is fundamental to many paleobiological studies (Simpson Reference Simpson1944; Stanley Reference Stanley1998; Levinton Reference Levinton2001). Although the significance of budding, anagenesis, and direct ancestry has been fiercely debated (Gould and Eldredge Reference Gould and Eldredge1977; Gould Reference Gould1980; Levinton and Chris Reference Levinton and Chris1980; Gingerich Reference Gingerich1985; Futuyma Reference Futuyma1987), these patterns have all been observed in the fossil record and can greatly impact inferences of evolutionary processes (MacLeod Reference MacLeod1991; Soul and Friedman Reference Soul and Friedman2017). This may be particularly important when constructing and testing hypotheses of phylogenetic relationships, in which accommodation of direct ancestry, whether as budding cladogenesis or anagenesis, may improve accuracy and yield otherwise unattainable insights (Gingerich Reference Gingerich, Cracraft and Eldredge1979; Fox et al. Reference Fox, Fisher and Leighton1999; Aze et al. Reference Aze, Ezard, Purvis, Coxall, Stewart, Wade and Pearson2011, Reference Aze, Ezard, Purvis, Coxall, Stewart, Wade and Pearson2013; Strotz and Allen Reference Strotz and Allen2013).
These issues are particularly relevant in taxa such as hominins (humans and all other taxa more closely related to humans than to chimpanzees), for which hypotheses of direct ancestry are often constructed and entertained qualitatively. Discerning human evolutionary patterns is a focus traceable to Darwin (Reference Darwin1871), and direct ancestry is frequently proposed (even if informally) for clusters of hominin fossils. Nonetheless, there are few formal routes for recognizing this status based on statistical evaluations of morphological traits.
Despite the importance of directly ancestral relationships in the fossil record, their role in phylogenetic reconstruction has remained underexplored. Stratocladistic methods were developed in part to explore use of temporal occurrence data to direct ancestor–descendant relationships (Fisher Reference Fisher2008). As a nonparametric approach, stratocladistics uses the criterion of maximum parsimony (MP) to minimize both the number of homoplasious character changes and unsampled stratigraphic intervals implied by topologies that include both bifurcating and serially linked segments of phyletic lineages. Several authors have expressed objections both to the use of temporal data in phylogenetic inference and the capability of available data to adequately test ancestral and collateral relationships between species (Smith Reference Smith2000). Nevertheless, direct ancestors are expected to occur in the fossil record (Foote Reference Foote1996), and the integration of temporal data has been shown to improve reconstruction accuracy over that of morphological analyses alone (Fox et al. Reference Fox, Fisher and Leighton1999). In addition, those earlier criticisms, which have been frequently raised by proponents of cladistic methodologies, are also less relevant when placed within the context of modern probabilistic approaches. While cladistics operates at the level of cladograms, both stratocladistics and recent probabilistic approaches reconstruct phylogenetic trees (Fisher Reference Fisher, Grande and Rieppel1994, Reference Fisher2008). This renders earlier criticism of temporal data largely irrelevant given the current landscape of phylogenetic methodology. As such, it stands to reason that probabilistic approaches are remiss when they fail to accommodate the possibility of ancestor–descendant relationships between taxa occurring at different times.
Two early attempts to explicitly extend the intent and logic of stratocladistics into parametric methods using maximum likelihood (ML) (Huelsenbeck and Rannala Reference Huelsenbeck and Rannala1997; Wagner Reference Wagner1998) showed strong potential, but did not fulfill all of the goals of stratocladistics, such as identifying direct ancestors. A more recent set of methods combine morphological, fossil preservation, and lineage diversification models to infer relationships and lineage divergence times in a Bayesian context (Pyron Reference Pyron2011). This framework has been extended to explicitly accommodate ancestor–descendant relationships by modeling lineage diversification and fossil preservation processes (Stadler Reference Stadler2010; Bapst and Hopkins Reference Bapst and Hopkins2017; Stadler et al. Reference Stadler, Gavryushkina, Warnock, Drummond and Heath2018), particularly through use of the “fossilized birth–death” (FBD) process (Heath et al. Reference Heath, Huelsenbeck and Stadler2014; Zhang et al. Reference Zhang, Stadler, Klopfstein, Heath and Ronquist2016). The FBD process models lineage diversification using speciation, extinction, and fossil preservation parameters, and therefore is very similar to previous “birth–death sampling” (BDS) models (e.g., Foote Reference Foote1997), but is usually estimated in a Bayesian context as an informative prior.
Although these approaches have been shown to be useful, several outstanding questions remain regarding their application to the fossil record. These methods have been largely developed for use in epidemiological systems, in which they are used to model molecular sequence evolution along single lineages. Although these systems are useful models of patterns in the fossil record, their sampling is often incomplete, and so they may sometimes call for different approaches in practice. Several a posteriori time-scaling (APT) approaches that accommodate ancestor–descendant relationships, such as cal3, also exist (Bapst Reference Bapst2013). These apply divergence times to unscaled cladograms that have been inferred from character data alone using a model similar to the FBD process. Although both Bayesian and APT approaches have been shown to be effective when applied to fossil taxa (Bapst Reference Bapst2013; Bapst and Hopkins Reference Bapst and Hopkins2017; Gavryushkina et al. Reference Gavryushkina, Heath, Ksepka, Stadler, Welch and Drummond2017; Stadler et al. Reference Stadler, Gavryushkina, Warnock, Drummond and Heath2018), there has been limited discussion of the statistical properties and identifiability of model parameters such as speciation and extinction rates, given incomplete fossil sampling. One more straightforward question regards the extent to which morphological data alone can provide evidence for direct ancestor–descendant relationships without modeling abstract lineage diversification parameters such as speciation or extinction rates. This may be especially important in clades with gap-prone fossil records, such as many lineages of terrestrial vertebrates. These taxa may provide less information from which to infer speciation and extinction parameters, and so a characterization of the evidence for hypotheses of direct ancestry provided by morphology alone will be important in the continued development of parametric approaches for phylogenetic inference in fossil species.
In this study, we describe an approach to phylogenetic inference that combines models of stratigraphic preservation and morphological evolution to reconstruct time-scaled phylogenies and distinguish between anagenesis and cladogenesis using ML and the Akaike information criterion (AIC). This approach seeks to simplify existing methods, such as the APT and FBD approaches described earlier, to clearly identify the signal for directly ancestral relationships presented by morphological data alone. Unlike more complex approaches, ours seeks only to identify information in morphological and temporal data that establishes differential support for ancestor–descendant and bifurcating relationships. As a result, our approach explicitly relies on models of morphological evolution and stratigraphic preservation rather than on models of lineage diversification.
We demonstrate the utility of our approach using the hominin fossil record. Hominins are a particularly compelling taxon for a case study. Although cladistic methods have been applied to hominin evolution (e.g., Delson et al. Reference Delson, Eldredge and Tattersall1977; Chamberlain and Wood Reference Chamberlain and Wood1987; Strait et al. Reference Strait, Grine and Moniz1997; Irish et al. Reference Irish, Guatelli-Steinberg, Legge, de Ruiter and Berger2013), prior authors’ exclusive focus on bifurcating relationships has precluded the study of direct ancestry, except in cases in which character polarity and specimen sampling have been carefully considered (Kimbel et al. Reference Kimbel, Lockwood, Ward, Leakey, Rak and Johanson2006). This limitation has been a major detriment to paleoanthropological studies, in which direct ancestry has often been considered to be especially important. The hominin fossil record has more recently been examined in a parametric context (Dembo et al. Reference Dembo, Matzke, Mooers and Collard2015, Reference Dembo, Radovčić, Garvin, Laird, Schroeder, Scott, Brophy, Ackermann, Musiba, de Ruiter, Mooers and Collard2016). However, these results still contain ambiguities and do not consider the possibility of direct ancestry. As a result, phylogenetic results have often conflicted with qualitative interpretations of hominin relationships. For instance, this divide occurs in previous treatments of Homo heidelbergensis and Homo antecessor. Previous studies have alternatively suggested that either H. heidelbergensis or H. antecessor are directly ancestral to Homo sapiens and Homo neanderthalensis (Mounier et al. Reference Mounier, Marchal and Condemi2009), while others have disagreed that H. antecessor is ancestral to H. heidelbergensis (Stringer Reference Stringer2012). Others have suggested that H. heidelbergensis is either a chronospecies directly ancestral to H. neanderthalensis or a direct ancestor shared by both H. sapiens and H. neanderthalensis (Rightmire Reference Rightmire1998; Rosas and Bermúdez De Castro Reference Rosas and Bermúdez De Castro1998; Stringer Reference Stringer2012). This uncertainty is underscored by the suggestion that cladistic methods and data are unreliable in their ability to describe hominin evolutionary patterns and history (Collard and Wood Reference Collard and Wood2000).
For our empirical exploration, we borrowed a morphological supermatrix from the literature (Dembo et al. Reference Dembo, Matzke, Mooers and Collard2015, Reference Dembo, Radovčić, Garvin, Laird, Schroeder, Scott, Brophy, Ackermann, Musiba, de Ruiter, Mooers and Collard2016). Our analysis identified several areas of direct ancestry in hominins, demonstrating cases in which temporal data both corroborate and refute results achieved using morphology alone. We recognize that increasingly thorough compilation of morphological data and new fossil discoveries is likely to continue to revise and refine current understanding of hominin evolution. Nevertheless, our approach sheds light on the capabilities of existing models and data to accommodate anagenesis and budding cladogenesis without relying on abstract diversification parameters. Moving forward, application of our method in hominins provides general demonstration of the importance of accommodating directly ancestral relationships in phylogenetic methods to generate deeper understanding of evolutionary patterns in the fossil record.
Materials and Methods
Use of Terms
We use the term “ancestor” modeled loosely after Gingerich (Reference Gingerich, Cracraft and Eldredge1979). Ancestors identified through our method represent taxa that display character states that are not sufficiently differentiated from those of a taxon occurring later in time to warrant assignment to a separate lineage. Operationally, we consider anagenesis as any evolutionary change occurring along these serially linked phyletic lineage segments. Strictly speaking, these instances might alternatively represent budding cladogenesis if new fossil evidence shows greater temporal overlap between ancestral and descendant taxa than is currently known. However, the gap-prone nature of much of the hominin fossil record complicates the ability to develop a precise mechanistic view of speciation patterns. Because candidate direct ancestors typically do not overlap in their temporal ranges with close relatives, it was simpler to consider all direct ancestors as connected to their descendants through “anagenesis” (again, given the lack of evidence for budding patterns in the stratigraphic record). Nevertheless, it would be straightforward to apply our approach to clades in which putative ancestors overlap with putative descendants.
Inference of ML Topology
Our approach evaluates the likelihood of candidate topologies using probabilistic models of fossil preservation and morphological evolution (Huelsenbeck and Rannala Reference Huelsenbeck and Rannala1997; Lewis Reference Lewis2001). We perform a semi-automated tree search by calculating the likelihoods of these models on a set of candidate topologies. This approach tests hypotheses of direct ancestry by combining branches with non-overlapping ranges and comparing cladogenetic and anagenetic models using the AIC. All code developed for these analyses is publicly available and implemented in the mandos package (www.github.com/carolinetomo/mandos).
Identification of a Fully Bifurcating ML Tree
Our approach has yet to be implemented with a fully automated tree-searching algorithm, so we combine semi-automated rearrangements with manual perturbations to search for the ML topology. This is done by exploring tree-space surrounding a starting tree estimated from morphological characters alone. In this study, we obtained an ML starting tree using RAxML v. 8.2.11, using the Mk model of morphological evolution (Lewis Reference Lewis2001). The morphological data were separated into partitions according to the number of possible states (i.e., binary, trinary, etc.) and analyzed under separate models. This partitioning scheme was maintained for all subsequent morphological likelihood calculations, including those used in the AIC comparisons presented in this study. We then performed a series of nearest neighbor interchange and subtree pruning and regrafting operations. These yielded a set of 1700 candidate topologies from which we identified the fully bifurcating topology best supported by both morphological and stratigraphic data. This tree provided a starting point drawn from both morphological and stratigraphic lines of evidence. From here, we explored ancestor–descendant relationships using the model-testing approach described in this paper.
Modeling Stratigraphic Preservation
Stratigraphic likelihoods were calculated under a homogeneous Poisson process of geologic preservation. When applied to phylogenetics, inference under this model has been shown to accurately recover simulated phylogenetic relationships (Huelsenbeck and Rannala Reference Huelsenbeck and Rannala1997). The likelihood function is derived from a Poisson process as the probability of observing the first occurrence (o f), last occurrence (o l) and number (n o) of occurrences in the stratigraphic record given some origination and extinction time (t f and t l, respectively) and preservation rate (λ). These likelihoods are calculated independently for each (i) of b total lineages and multiplied to yield the overall tree likelihood (see Huelsenbeck and Rannala [1997] for full derivation):
$$\hat{L} = \mathop \prod \limits_{i = 1}^b \displaystyle{{{\lpar {o_f^i -o_l^i} \rpar }^{n_o^i -2} {\rm \lambda} ^{n_o^i} e^{-\lambda \lpar {t_f^i -t_l^i} \rpar }} \over {(n_o^i -2)!}}$$This equation reaches its maximum as branching and extinction times approach the first and last occurrences in the fossil record, and so the likelihood is maximized across the tree when the total amount of unsampled time implied by the topology is minimized. Although this approach differs from stratocladistic parsimony in its treatment of occurrence data as continuous rather than discrete, this property causes the preservation model to behave as a statistical formalization of the stratigraphic parsimony debt calculations undertaken in parsimony-based stratocladistic analyses. We estimate lineage origination and termination times along with preservation rate using multivariate numerical optimization routines implemented in SciPy (Jones et al. Reference Jones, Oliphant and Peterson2001). When combined with morphological data and models of character evolution, this approach represents a comprehensive extension of traditional stratocladistics based entirely on probabilistic models.
Identification of Direct Ancestry
Using the bifurcating topology with the highest likelihood under both stratigraphic and morphological models, we manually identified a set of potential ancestor–descendant relationships. To explore a more comprehensive range of both anagenetic and cladogenetic arrangements, we also extensively perturbed results manually and compared likelihoods. Although a fully automated tree-searching approach will ultimately be desirable in future versions of our method, our approach to tree-searching is similar to those used in previous stratocladistic studies. Starting with the fully bifurcating tree, we identified putative ancestor–descendant arrangements by collapsing each branch with a temporal range beginning earlier than the range represented by its sister lineage. We isolated each putative episode of anagenesis and compared the morphological and stratigraphic likelihoods of anagenetic and cladogenetic arrangements using AIC scores. This was required, because cladogenetic nodes assume one more parameter than anagenetic nodes (i.e., the branch length or node height connecting the new lineage), and so bifurcating trees contain more parameters than anagenetic trees. Comparison of AIC scores enables a comparison of the relative quality of models with different numbers of parameters. AIC score is calculated from the number of model parameters (k) and the likelihood (L):
$${\rm AIC}\, = \,2k-2\ln L$$More complex models almost always have higher likelihoods than simpler models, because added parameters allow a better fit to slight deviations in the data. Because they represent a more parameter-rich phylogenetic model, likelihoods from fully bifurcating trees are not directly comparable to ancestor–descendant trees. AIC scores represent the amount of information lost by a model when representing data, with lower scores being indicative of models that preserve a greater amount of information. AIC accommodates differences in parameter count between models by penalizing the addition of new parameters, seeking to optimize the trade-off between improvements in model fit associated with added parameters and the loss of statistical power that results from overparameterization. In our use, AIC facilitates comparison of phylogenetic models of differing dimensions by penalizing the addition of branches that are better explained through an anagenetic pattern. The parameter count k for each phylogeny is calculated by summing the number of estimated branch lengths for each tree with the number of parameters used in the partitioned Mk substitution model.
Under the Poisson preservation model, stratigraphic likelihood predictably improves when unsampled time implied by the phylogeny is reduced, and so the acceptance of ancestor–descendant arrangements also requires the support of morphology. We accommodated direct ancestry using a novel calculation of morphological likelihood wherein the probability of transitioning from an observed, rather than an uncertain, parental character state to a single or multiple descendant character states is calculated under the Mk model (Fig. 1). This calculation differs from that used on multifurcating nodes. Because the sequences at internal nodes representing unobserved taxa are unknown, the likelihood of character data at the tips is typically calculated by summing over all possible states at each unobserved internal node (Felsenstein Reference Felsenstein1981). However, when dealing with direct ancestry, the sets of character states at some internal nodes are known and the marginal likelihoods for states absent in the ancestor are necessarily zero. In these cases, the likelihood is calculated as the conditional probability of observing the set of traits at the tips given the set of traits possessed by the putative ancestor. Ancestor–descendant arrangements are only accepted when there has been a sufficiently small amount of character change. This minimum threshold varies depending on the informativeness of the data set and the trade-off between model fit and parsimony that are weighted when using AIC. This procedure resembles model-testing procedures used to reconstruct ancestral DNA sequences, which compare conditional likelihoods of different permutations of character states at ancestral nodes (Yang et al. Reference Yang, Kumar and Nei1995). A similar effect could also be achieved by fixing the branch length leading to the putative ancestor to zero and performing the standard pruning algorithm (Felsenstein Reference Felsenstein1981). Although the two calculations are equivalent, our approach treats direct ancestry explicitly in the likelihood calculation and tree representation. Our use of trait evolution models enables morphology to occupy a central role in identifying direct ancestors. As a result, this approach can, in principle, be applied without temporal data (although we do not recommend this practice). A demonstration and test of this method using simulated data are provided in the Supplementary Material (https://github.com/carolinetomo/hominin_anagenesis). These preliminary validations of the implementation show that the approach identifies direct ancestors correctly, even when rates of change are very uneven across taxa. However, the approach can incorrectly prefer arrangements in which a taxon with a low rate of morphological change relative to its sister taxon is collapsed into a directly ancestral position. This result demonstrates the importance of surveying broadly across the anatomy when constructing character matrices and the key role played by temporal data in constraining a set of possible ancestor–descendant relationships.
Comparison between anagenetic and cladogenetic trees. A, Likelihoods are calculated on bifurcating arrangements using the standard “pruning algorithm.” B, Direct ancestor likelihood calculation.
Morphological Matrix
We performed our analysis on a supermatrix of 391 discrete craniodental characters compiled by Dembo and colleagues (Dembo et al. Reference Dembo, Matzke, Mooers and Collard2015, Reference Dembo, Radovčić, Garvin, Laird, Schroeder, Scott, Brophy, Ackermann, Musiba, de Ruiter, Mooers and Collard2016). We removed all ambiguous character states, as the researchers did not identify whether these were truly ambiguous or polymorphic. While ambiguous character codings are unlikely to provide significant phylogenetic information, existing Markov models of discrete character evolution do not accommodate polymorphism. We excluded the taxa Kenyanthropus platyops and Homo naledi from the present analysis. The features that are diagnostic of K. playtyops have been suggested to result from taphonomic distortion resulting from matrix expansion, rather than from true biological differences (White Reference White2003). Thus, we omitted this taxon in hopes of shedding greater light on the remaining, more widely accepted hominin taxa. Homo naledi was omitted because the data provided in the original study yielded an ML topology placing H. naledi as sister to H. sapiens. Although the phylogenetic affinity of H. naledi is a major outstanding question in paleoanthropology, the confusing signal presented by the H. naledi data, which are relatively recently acquired and therefore represent less well-studied fossils overall, reduced our confidence in the ability of this data set to resolve its placement. Therefore, to avoid any confounding effects from the potential unreliability of the H. naledi data, we performed our analyses on the remaining subset of the data after H. naledi was removed. This enabled us to explore the phylogenetic relationships between better-known hominins.
Geologic Occurrence Times
We surveyed the literature to obtain the observed temporal range of each taxon in continuous time. Reported radiometric dates for the oldest and youngest fossils were taken as the first and last observations. Some specimens are ambiguous in their taxonomic assignment; these were excluded from the analysis. We also gathered the number of total occurrences as the number of localities where each taxon has been identified as listed in MacLatchy et al. (Reference MacLatchy, Desilva, Sanders, Wood, Sanders and Werdelin2010) and supplemented these with additional localities identified in the literature. Cases in which multiple specimens belonging to the same taxon have been identified at a single locality were treated as single occurrences. Although we recognize the potential ambiguity in delineating between sites, localities, and occurrences, we attempted to coarsely characterize the total number of occurrences using the number of sites at which each taxon occurs. This approach is more likely to underestimate the number of occurrences than overestimate them, which we expect to yield more conservative statistical support for competing topologies under the preservation model. A comprehensive list of the sites used to define temporal ranges for all taxa is provided in the Supplementary Material.
Results and Discussion
Anagenesis in the Hominin Fossil Record
Our analysis yielded evidence for several instances of anagenesis in the hominin fossil record (Fig. 2). It reconstructed Australopithecus anamensis as directly ancestral to Australopithecus afarensis. This result agrees with broad acceptance of Au. anamensis and Au. afarensis as phyletically linked chronospecies (Leakey et al. Reference Leakey, Feibel, McDougall and Walker1995; Ward et al. Reference Ward, Leakey and Walker2001; Kimbel et al. Reference Kimbel, Lockwood, Ward, Leakey, Rak and Johanson2006). Although our analysis recovered Ardipithecus ramidus as sister to the anagenetic Au. anamesis-Au. afarensis branch, the morphological data did not support the collapse of Ar. ramidus. Nevertheless, Ar. ramidus possessed poor character sampling in the matrix, and so its placement should be regarded as tentative. Sahelanthropus tchadensis is recovered as a direct ancestor to the rest of the hominin clade. This result should also be treated cautiously due to the small number of characters recovered for analysis of this portion of the phylogeny, but it is in line with Sahelanthropus’s status as the oldest recognized hominin (Brunet et al. Reference Brunet, Guy, Pilbeam, Mackaye, Likius, Ahounta, Beauvilain, Blondel, Bocherens, Boisserie, De Bonis, Coppens, Dejax, Denys, Duringer, Eisenmann, Fanone, Fronty, Geraads, Lehmann, Lihoreau, Louchart, Mahamat, Merceron, Mouchelin, Otero, Campomanes, De Leon, Rage, Sapanet, Schuster, Sudre, Tassy, Valentin, Vignaud, Viriot, Zazzo and Zollikofer2002; Guy et al. Reference Guy, Lieberman, Pilbeam, de Leon, Likius, Mackaye, Vignaud, Zollikofer and Brunet2005; Zollikofer et al. Reference Zollikofer, Ponce de León, Lieberman, Guy, Pilbeam, Likius, Mackaye, Vignaud and Brunet2005).
Phylogenetic relationships between hominin species. A, Reconstruction with the best Akaike information criterion (AIC) score when temporal data are considered alongside morphological data and anagenesis is accommodated. B, C, Areas of the tree inferred from morphology alone that differ from A. Both regions display misleading results when anagenesis is not considered. This is reflected in the improvement in AIC observed in the preferred anagenetic topology. Silhouettes obtained from phylopic.org.
Consistent with an early appraisal (Asfaw et al. Reference Asfaw, White, Lovejoy, Latimer, Simpson and Suwa1999), our final analysis inferred Australopithecus garhi to be directly ancestral to the Homo clade (Fig. 2A). This conflicts with cladistic analyses that placed Au. garhi as an outgroup to Au. africanus, Paranthropus, and Homo (Strait Reference Strait1999). However, when anagenesis is not considered and phylogeny is inferred from morphology alone, we recover the same placement for Au. garhi as the cladistic result (Fig. 2B). This may reflect the constraint that strictly bifurcating methods impose on phylogenetic reconstructions among fossil taxa. However, preference for Au. garhi as ancestral to Homo is relatively weak, with an AIC score (9264.3) that is only slightly better than that of the Au. garhi outgroup hypothesis (9271.11).
Our results also better reconcile quantitative and qualitative interpretations of the evolutionary relationships (involving potentially some combination of direct ancestry and collateral connections) among the Homo species leading to modern humans and Neanderthals. Paleoanthropologists have variously interpreted Homo heidelbergensis as either: (1) a direct ancestor to both Neanderthals and modern humans (Rightmire Reference Rightmire1998; Mounier et al. Reference Mounier, Marchal and Condemi2009; Stringer Reference Stringer2012; Buck and Stringer Reference Buck and Stringer2014) or (2) a chronospecies leading to Neanderthals (Rosas and Bermúdez De Castro Reference Rosas and Bermúdez De Castro1998). However, quantitative phylogenetic analysis has placed H. heidelbergensis as sister to Neanderthals (Dembo et al. Reference Dembo, Matzke, Mooers and Collard2015). Our initial result based on morphological characters alone places H. heidelbergensis as sister to H. sapiens. However, AIC support improves substantially when H. heidelbergensis is collapsed to represent a direct ancestor that preceded the split between modern humans and Neanderthals (Fig. 3). We also uncover one other instance of anagenesis in this clade. Opinions are divided as to whether H. antecessor represents a direct ancestor of later hominin species or is an evolutionary dead end (Bermúdez de Castro et al. Reference Bermúdez de Castro, Arsuaga, Carbonell, Rosas, Martínez and Mosquera1997; Stringer Reference Stringer2012; Dembo et al. Reference Dembo, Matzke, Mooers and Collard2015), but our results provide support for combining H. antecessor and H. heidelbergensis into a single lineage. This suggests a long episode of anagenetic evolution immediately before the divergence between modern humans and Neanderthals.
A–F, Akaike information criterion (AIC) scores calculated for each possible arrangement between Homo sapiens, Homo neanderthalensis, and Homo heidelbergensis.
Overall, our reconstruction of relationships among later species of Homo immediately preceding and encompassing modern humans and Neanderthals supports the hypothesis that H. heidelbergensis is directly ancestral to both modern humans and Neanderthals. This result differs from the hypothesis supported by the analysis performed by Dembo and colleagues (Reference Dembo, Matzke, Mooers and Collard2015), and is instead more consistent with an earlier exploratory statistical analysis (Mounier et al. Reference Mounier, Marchal and Condemi2009) and with the position frequently suggested by paleoanthropologists (Rightmire Reference Rightmire1998; Stringer Reference Stringer2012; Buck and Stringer Reference Buck and Stringer2014). While previous phylogenetic analyses have yielded results that equivocate or disagree with the common interpretation of H. heidelbergensis as the last common ancestor of modern humans and Neanderthals, our analysis shows that the consideration of direct ancestry can generate statistical support for phylogenetic results that conform more closely to positions generated through researchers’ subjective interpretations based on morphological considerations and/or stratigraphic and biogeographic assessments and on exploratory data analyses. This finding supports a general argument against the use of cladistic and phylogenetic methods that are restricted to bifurcating relationships in fossil taxa, in which the possibility of variability in evolutionary mode (i.e., occurrence of both anagenesis and cladogenesis) is at odds with an assumption that evolution proceeds by cladogenesis alone (Fig. 2).
Our results differ markedly from previous phylogenetic studies seeking to reconstruct hominin phylogeny using probabilistic and cladistic methods. In key regions of the tree, results achieved under our method reveal support for hypotheses more consistent with many qualitative interpretations of hominin relationships, demonstrating the importance of explicitly accommodating anagenesis in the phylogenetic reconstruction of fossil taxa. This may explain some of the historical difficulty in reconciling paleontological interpretations of hominin relationships with cladistic results. For instance, Dembo and colleagues’ (Reference Dembo, Matzke, Mooers and Collard2015) results are inconsistent with earlier suggestions that Au. anamensis and Au. afarensis are chronospecies differentiated through anagenesis. However, by considering ancestor–descendant relationships and incorporating temporal data, our analysis reveals that a linked Au. anamensis-Au. afarensis lineage is the arrangement most strongly supported by the data. Generally, this suggests that cladograms and strictly bifurcating phylogenies may be inadequate when describing evolutionary relationships between fossil taxa, and so ancestor–descendant relationships should be considered during topological inference. Further, we argue that explicit testing of ancestor–descendant relationships is important, even in cases in which bifurcating trees are not wholly misleading, as their omission precludes us from considering the full range of possible evolutionary scenarios. Although previous studies addressing hominin phylogeny using probabilistic methods represent a significant step forward in weighing alternative evolutionary hypotheses, we suggest that their phylogenetic reconstructions have suffered from methodological limitations that were not generally perceived. We suggest that overcoming these limitations can provide a substantial step forward in closing the gap between the paleobiologists’ interpretations and previous cladistic and phylogenetic results. In particular, we show through our analyses that the apparent discordance between quantitative and qualitative assessments of evolutionary relationships can be reconciled by extending phylogenetic models to explicitly accommodate anagenesis.
Ancestors, Anagenesis, and Evolutionary Processes
Our method does not seek to distinguish between speciation modes at a mechanistic level. As noted by Fisher (Reference Fisher2008), overlap between the extinction and origination times of taxon pairs does not necessarily preclude an ancestor–descendant relationship. Instead, it is possible for ancestral and descendant lineages to coexist. This reality complicates the identification and evolutionary interpretation of ancestor–descendant relationships from temporal data alone when the observed ranges of putative ancestor–descendant pairs do not overlap. Detailed interpretations would require distinction between anagenesis and budding. This requires BDS models, which contain additional lineage diversification (origination and extinction) parameters. Such models are currently implemented in the cal3 time-scaling method, which seeks to distinguish among splitting, budding, and anagenesis from a cladogram with a fixed topology (Bapst and Hopkins Reference Bapst and Hopkins2017). Our approach is distinct both theoretically and operationally from cal3, instead relying most heavily on morphological evidence to weigh the likelihood of ancestor–descendant relationships without considering the completeness of fossil sampling. At the population level, one might expect some temporal overlap between ancestral and descendant taxa undergoing anagenetic change, especially if the participants are widely distributed geographically (for additional remarks on the relevance of biogeographic data see Fisher Reference Fisher, Grande and Rieppel1994: pp. 158–164; Reference Fisher2008: p. 380). As a result, even if the temporal ranges corresponding to taxa identified as ancestor–descendant pairs are discovered to underestimate slightly the time of coexistence, the inference may simply highlight the fuzziness in the taxonomic placement of fossils belonging to lineages undergoing continuous transformation and in discerning between anagenesis and evolutionary budding given incomplete sampling (Fig. 4). This interpretation is consistent with previous authors’ treatment of temporal ranges when identifying anagenesis between taxa, which has allowed a period of overlap between putative ancestor–descendant pairs (Aze et al. Reference Aze, Ezard, Purvis, Coxall, Stewart, Wade and Pearson2011; Strotz and Allen Reference Strotz and Allen2013).
A, Speciation mode interpreted as budding when sampling is complete. B, Incomplete stratigraphic sampling may create an inability to distinguish between anagenesis and budding when sampling is sparse.
Under our method, ancestor–descendant relationships might represent either true anagenesis (i.e., a single population undergoing gradual transformation) or some form of budding cladogenesis. Previous researchers have argued that true anagenesis is rare compared with budding when analyzing the fossil records of densely sampled marine invertebrate lineages using more complex preservation models (Bapst and Hopkins Reference Bapst and Hopkins2017). Nevertheless, we suggest that distinction between these two modes may often be impossible in terrestrial vertebrate lineages with large sampling gaps. For example, our results among early hominin species include multiple inferred direct ancestors, but the large gaps in stratigraphic sampling throughout this region of the tree hamper the ability to determine whether these relationships represent true anagenesis or budding that has been obfuscated by poor sampling.
As employed here, our method makes no attempt to distinguish between budding and anagenesis. Although temporal data occupy an important place in our approach, they are largely used as a guide to constrain the set of possible ancestors and descendants and to provide additional insight when morphological data are equivocal. Coarsely speaking, our approach focuses on patterns in morphological differentiation and lineage disparification, while approaches such as cal3 model lineage diversification. FBD models employed in Bayesian analyses weigh both morphological data and diversification parameter estimates and so combine elements of both approaches. Further empirical and simulation-based work is needed to determine the differences in behavior between methods that test ancestor–descendant relationships in a BDS/FBD framework, and ours. We speculate that their relative accuracy may depend largely on the completeness of sampling in the rock record and the correlation strength between morphological change and lineage diversification, although other factors may also be important.
The scales at which phylogenetic data are sampled may further complicate mechanistic evolutionary interpretations. Morphological character matrices often lack samples across the entire stratigraphic range of each taxon, so in the absence of evidence to the contrary, analysts often assume morphological stasis within lineages. It is therefore often impossible to observe gradual morphological change within and between taxa. These considerations might cause anagenetic relationships identified here to represent either true anagenesis or some form of “pseudo-anagenetic” budding in which stratigraphic and morphological data appear consistent with anagenesis but the persistence of the ancestor has not been sampled. The ancestor–descendant relationships identified by our method may be interpreted in several ways. As sampled, these results may be roughly conceived as anagenesis in the sense that the mode of evolutionary change between taxa is indistinguishable from evolution occurring along a single lineage, depending upon the completeness of sampling and the degree to which morphological disparity correlates with true biological or evolutionary species diversity. This interpretation is consistent with historical usage by paleobiologists (Gingerich Reference Gingerich, Cracraft and Eldredge1979; Smith Reference Smith1994; Levinton Reference Levinton2001). Thus, our approach seeks to reveal the extent to which existing cladistic and temporal data can provide evidence for nonbranching evolutionary modes and does not seek to resolve conceptual issues that may stem from incomplete sampling, lineage diversification, or population-level evolutionary change. Regardless of the fine-scale evolutionary interpretation, failure to accommodate phyletic change and ancestor–descendant relationships when inferring phylogenetic relationships can generate views of evolutionary history that are positively misleading, in the sense that inaccurate results do not improve with added data.
Some Practical Methodological Considerations
Concerns regarding the accuracy of probabilistic approaches have been raised, stemming from the reliance of these methods on the overly simplistic Lewis Mk model of morphological evolution (Goloboff et al. Reference Goloboff, Torres and Arias2017). These critics advocate the use of cladistic methods, arguing that Markov models inadequately capture the complexities of morphological evolution. Although we agree that existing substitution models oversimplify these processes, our results suggest that the accommodation of ancestor–descendant hypotheses in probabilistic methods can improve the fidelity of phylogenetic reconstructions, even when Lewis Mk is used as the underlying model of morphological change. As a result, concerns regarding the adequacy of existing morphological substitution models may be partially alleviated by considering hypotheses of direct ancestry. This is supported by simulation work showing that stratocladistics outperforms cladistics in topology reconstruction (Fox et al. Reference Fox, Fisher and Leighton1999). Further exploration is needed to demonstrate more thoroughly the limitations of our new approach, which builds upon stratocladistics by incorporating the benefits of probabilistic analyses, including (1) more explicit statements of the assumptions involved and (2) the ability to weigh competing models using modern inferential criteria.
The method we describe seeks to enhance understanding of the fossil record by explicitly testing support for existing hypotheses of direct ancestry while attempting to make simpler assumptions than stratocladistics or recently developed Bayesian methods. Although Bayesian implementations that employ the FBD model also seek to identify direct ancestry (Zhang et al. Reference Zhang, Stadler, Klopfstein, Heath and Ronquist2016), they often struggle to parse through complex signals to choose between competing ancestor–descendant and collateral hypotheses (Luo et al. Reference Luo, Duchene, Zhang, Zhu and Ho2018). Although future extensions to those approaches may improve their response in this context, we view our approach as a foundational, minimally complex framework for exploring the behavior of probabilistic models when evaluating support for direct ancestry in temporal and morphological data. As such, our method should be viewed as a complement to, rather than a simplification of, existing Bayesian approaches. Our method encourages examination of the informativeness of the data without the increased complexity of assessing prior probabilities conditioned upon models of lineage diversification. Thus, our method differs from both existing Bayesian and parametric APT approaches by explicitly omitting diversification parameters and instead placing morphological data in a central role when evaluating hypotheses of direct ancestry. In our approach, temporal data help to delineate the set of possible ancestors and play an important role in measuring the fit of candidate trees to the observed stratigraphic record. Our method does not seek to reconstruct diversification processes, and instead focuses on identifying hypotheses that best describe only the information contained within morphological and temporal data sets. Assessment of information contained within data sets and tests of hypotheses can also be achieved using Bayesian approaches (Lewis et al. Reference Lewis, Chen, Kuo, Lewis, Fučíková, Neupane, Wang and Shi2016), but likelihoodist approaches such as ours streamline these procedures by reducing complications presented by prior probabilities.
Although Bayesian methods can be beneficial in certain circumstances, our method simplifies identification of anagenetic hypotheses using evolutionary and stratigraphic models. We observe that the likelihood surface surrounding certain nodes may possess low peaks (Fig. 2), which likely results from sparse sampling and relatively low information in the morphological characters. Because Bayesian approaches often average results across this surface, it is possible that they may fail to capture those relationships best supported by the data by including information from weakly supported hypotheses. This is of greater concern in paleontological than neontological data, because the increased abundance of molecular data is often likely to result in more clearly defined peaks in the likelihood surface. In cases in which information is sparse, likelihoodist approaches such as ours offer the benefit of filtering through noisy and equivocal signal to reveal the hypothesis best supported by the data. Although these benefits are also achievable through Bayesian approaches, additional caution must be taken to select priors that do not dominate weakly informative data. In addition, careful thought should be given when summarizing the posterior/likelihood surface. Averaging across a relatively flat surface might yield poor results (Yang and Zhu Reference Yang and Zhu2018), while the comparison of individual point estimates, as is done here, may more clearly shed light on best-supported models while clearly contrasting this support relative to competing models. Moving forward, further extensions to our method that more clearly evaluate uncertainty across the likelihood surface will be useful.
How Can We Proceed?
Anagenesis and splitting cladogenesis were most pertinent to our analysis of hominin evolution, and the straightforward dichotomy between these simplified the assumptions and interpretation of our tests. Our approach may need to explicitly accommodate budding before being applied to groups with very dense fossil records. However, the assumptions required, which may often include morphological stasis within lineages, may make application to some phylogenetic data sets impractical. This is especially true in cases such as hominins, for which the fossil record implies large sampling gaps, and characters representative across stratigraphic ranges are often sampled from only a single individual. Expansions of our approach through implementation of new models will further test the implications of existing paleontological data sets for reconstructing complex evolutionary and geologic processes over deep timescales. We hope that the example provided here will encourage integration of more diverse evolutionary modes into phylogenetic methods, yielding better explanations of temporal patterns in critical parts of life's history.
As we emphasized earlier, the approach described here should be viewed as an attempt to explore the capability of phylogenetic methods to identify the signature of direct ancestry using existing models to interrogate morphological and stratigraphic data. In doing so, we acknowledge that there are many complicated biological and geological factors that could be incorporated into this framework. For instance, previous researchers have accommodated heterogeneity in fossil preservation rates across time and among lineages (Foote Reference Foote1997, Reference Foote2001; Gavryushkina et al. Reference Gavryushkina, Welch, Stadler and Drummond2014). There have also been several concerns raised in the literature regarding the adequacy of existing models of discrete trait evolution to inform complex evolutionary scenarios (Brown et al. Reference Brown, Parins-Fukuchi, Stull, Vargas and Smith2017; Goloboff et al. Reference Goloboff, Torres and Arias2017). Alternative models that use continuous characters may help to improve some of these issues (Felsenstein Reference Felsenstein1988; Goloboff et al. Reference Goloboff, Mattoni and Quinteros2006; Parins-Fukuchi Reference Parins-Fukuchi2017). Moving forward, elaborations making use of new data sources and models will only continue to improve resolution of evolutionary patterns in the fossil record.
Finally, we acknowledge that our empirical results beg qualification. In particular, we expect that future studies will generate a more comprehensive and authoritative view of hominin evolution as improved data continue to become available. For instance, although we are currently cautious about making strong statements concerning the ancestral position of Sahelanthropus using this data set, additional information may resolve this issue. This may be the case for several other areas of the hominin tree, which may be better resolved as temporal and taxonomic gaps in sampling are better filled by new discoveries. In addition, we concede the possibility that more comprehensive automated tree-searching routines may reveal support for hypotheses that we failed to consider under our semi-automated approach. Therefore, instead of providing an authoritative view of hominin evolution, our study provides a springboard for future studies by showing that the accommodation of anagenesis can improve our view of the processes and relationships underpinning the evolution of fossil taxa. Nevertheless, due to the improved support and congruence of hypotheses that explicitly consider ancestor–descendant relationships, we recommend that future phylogenetic studies in hominins avoid methods that only consider bifurcating relationships. Future studies that build upon existing work in other taxa will also be important to better characterize the extent to which this suggestion can be generalized across the tree of life. Although the accommodation of directly ancestral relationships is especially relevant in hominin taxa, for which hypotheses of anagenesis have been long entertained through qualitative anatomical assessment, these results may also be important in other taxa. Further empirical work will be needed to develop a better understanding of the extent to which the consideration of direct ancestors can improve resolution of evolutionary patterns throughout the fossil record.
Acknowledgments
C.P.-F. would like to thank D. Bapst, J. W. Brown, C. W. Dick, M. Friedman, S. A. Smith, and G. W. Stull for helpful discussion of the article and method. We would also like to thank P. Wagner and one anonymous reviewer for comments that helped refine the manuscript.
Author Contributions
C.P.-F. designed and implemented the ML “stratophylogenetic” procedure, performed analyses, interpreted results, and drafted the article. D.C.F. made conceptual contributions to the approach introduced here, contributed to the article, and helped interpret results. E.G. gathered temporal occurrence ranges, helped interpret results, and edited the article. L.M.M. helped focus the scope of the hominin analysis and interpret results and contributed to the article.
Supplementary Data
Software used in this study is implemented in the mandos package and is available on GitHub (github.com/carolinetomo/mandos). All trees tested using the stratocladistic likelihood procedure and their associated AIC scores are available on GitHub (https://github.com/carolinetomo/hominin_anagenesis). This repository also contains the morphological matrix and stratigraphic range files and all scripts needed to reproduce results. Stratigraphic occurrences and their associated citations are provided as Supplementary Material. The Supplementary Material also includes an online appendix that describes and tests the morphological likelihood calculation on trees containing ancestor–descendant relationships. Data matrices and mandos scripts used to run analyses are also stored on Dryad (doi:10.5061/dryad.6d2h8bp).
