Continuous-Time versus Punctuated Change in Phylogenetic Models

Different possible macroevolutionary scenarios often predict different relationships between phylogeny and observable data. Models of character evolution exemplify this. The most common method used for anatomical data is Lewis’s (Reference Lewis2001) Mk model. Here, the likelihood of a specific combination of relationships, diversification times, and rates of change is the Poisson probability of net stasis or, ultimately, transitioning to another state is given by:

(1) $$ {e}^{Q_C} $$

where Qc in its simplest form is a rate matrix:

(2) $$ {Q}_C=t\times \left(\begin{array}{cccc}-\alpha & \frac{\alpha }{k-1}& \dots & \frac{\alpha }{k-1}\\ {}& & & \\ {}\frac{\alpha }{k-1}& -\alpha & \dots & \frac{\alpha }{k-1}\\ {}& & & \\ {}\vdots & \vdots & \ddots & \vdots \\ {}& & & \\ {}\frac{\alpha }{k-1}& \frac{\alpha }{k-1}& \dots & -\alpha \end{array}\right) $$

with t being the amount of time over which a character might change (e.g., the span between a divergence time and the first appearance of a taxon), α is the instantaneous (Poisson) rate of change, and αt therefore is the expected amount of change. (It is unnecessary to place t outside the matrix and using −αt or − $ \frac{-\alpha t}{k-1} $ in each cell will yield the same results; we use this form because phylogenetic analyses vary t and α independently and because the scripts for performing these analyses require separating α and t.) Following matrix exponentiation, the probability of net stasis is $ \frac{1+\left(k-1\right){e}^{- k\alpha t}}{k} $ , and the probability of any transition is $ \frac{1-{e}^{- k\alpha t}}{k} $ , where k is the number of character states (e.g., 3 in eq. 2; Lewis Reference Lewis2001). Note that the probability of net stasis is not the Poisson probability of no change (e ^−αt ) but instead is the probability of no net change: that is, starting at State X and ending at State X. This includes both no character change and 2+ characters changes ultimately leading to reversal back to the original condition. Thus, for a binary character, the probability of net stasis is the Poisson probability of zero changes plus the summed probabilities of all even numbers of change. Correspondingly, the probability of “ultimate” transition for a binary character is the summed probability of all odd numbers of changes. When there are 3+ states, these probabilities are weighted by the probability that (say) 3 changes will lead back to the original state or from State X to State Y. More complex models exist that allow for biased state transitions (including driven trends sensu McShea [Reference McShea1994]), variation in probabilities of state transitions over phylogeny (Nylander et al. Reference Nylander, Ronquist, Huelsenbeck and Nieves-Aldrey2004; Wright et al. Reference Wright, Lloyd and Hillis2016), and correlated character change (e.g., Billet and Bardin Reference Billet and Bardin2018). The Mk model itself is a special case of all of these more complex models in which state transition rates are equal under all circumstances (Guimarães Fabreti and Höhna Reference Guimarães Fabreti and Höhna2023). Thus, these variants all assume that anatomical change can happen at any point in time rather than being concentrated in speciation events.

Because punctuated models posit that change is limited to cladogenetic events, we might expect a punctuated analogue to the Mk model to use binomial/multinomial probability. However, this would require that we know the number of cladogenetic events that occur over some duration t. If cladogenesis itself is a continuous-time process with rate λ, then the expected number of events over time t is λt. If the probability of change per branching event is ε, then the expected amount of change over time t is simply ελt. This led Lewis (Reference Lewis2001: p. 917) to suggest that the Mk model accommodates punctuated change, making no assumptions about how the change is distributed across a given branch. Because this ignores information about necessary cladogenetic events that result in clades having 2+ contemporaneous lineages, ελt predicts punctuated anagenetic change in which change is neither continuous over time nor attached to cladogenetic events (Wagner and Marcot Reference Wagner, Marcot, Alroy and Hunt2010). This differs from Eldredge and Gould’s model, which posits that cladogenesis and anatomical change have a common cause. Phylogenetic topologies will imply some minimum number of cladogenetic events, m, which means that expected change is ε(m + λt). We can assume m = 0 only for a possible anagenetic topology (Fig. 1A). Given a budding cladogenetic topology (Fig. 1B), m = 1 for the branch leading to the daughter species, and expected change is ε(1 + λt). Given a bifurcating pattern between two true sister taxa, there are two possibilities for any one branch (Fig. 1C,D). There must be one branching event, but we usually will have no reason to think that the “necessary” branching event leads to either branch: it might be either a continuation of the ancestral line or it might be the product of cladogenesis. Thus, the prior probability that the necessary cladogenetic event leads to the branch in question is 0.5, and expected change now is: $ \frac{\left(\varepsilon \left[1+\unicode{x03BB} t\right]\right)+\left(\varepsilon \left[0+\unicode{x03BB} t\right]\right)}{2}=\varepsilon \left(0.5+\unicode{x03BB} t\right) $ . This effectively makes the minimum number of branching events leading to either of two true sister taxa m = 0.5, and the expected number over branch duration t = 0.5 + λt. A key difference that emerges is that even if t = 0, the expected change is greater than zero given any phylogenetic topology in which there must have been at least one cladogenetic event (i.e., m = 0.5 or 1.0).

Figure 1. Four possible cladogenetic histories leading to taxon X. Thin black vertical lines denote necessary cladogenesis. Blue lines denote time over which continuous change and/or additional cladogenetic events might occur. A, X is the descendant of an observed species that last appears before X first appears. Here, there might be no cladogenetic events (m = 0). B, X is the descendant of a species that co-occurs with X. Here, there must be at least one cladogenetic events (m = 1). Bifurcation, where neither X nor its sister taxon are ancestral to one another. There must be one cladogenetic event; however, C, X might be the continuation of the original ancestor and thus need not have any cladogenetic events separating it from its common ancestor with X’s sister taxon. D, Alternatively, X represents the line derived from the cladogenetic event. Assuming both possibilities to be equiprobable, here m = 0.5.

Following Wagner and Marcot (Reference Wagner, Marcot, Alroy and Hunt2010), a “true” punctuated equilibrium analogue to the Mk model in its simplest form is given by:

(3) $$ {Q}_P=t\times \left(\begin{array}{cccc}-\varepsilon \left(\frac{m+\lambda t}{t}\right)& \varepsilon \left(\frac{m+\lambda t}{\left(k-1\right)t}\right)& \dots & \varepsilon \left(\frac{m+\lambda t}{\left(k-1\right)t}\right)\\ {}& & & \\ {}& & & \\ {}\varepsilon \left(\frac{m+\lambda t}{\left(k-1\right)t}\right)& -\varepsilon \left(\frac{m+\lambda t}{t}\right)& \dots & \varepsilon \left(\frac{m+\lambda t}{\left(k-1\right)t}\right)\\ {}& & & \\ {}\vdots & \vdots & \ddots & \vdots \\ {}& & & \\ {}\varepsilon \left(\frac{m+\lambda t}{\left(k-1\right)t}\right)& \varepsilon \left(\frac{m+\lambda t}{\left(k-1\right)t}\right)& \dots & -\varepsilon \left(\frac{m+\lambda t}{t}\right)\end{array}\right) $$

with the probability of net stasis being $ \frac{1+\left(k-1\right){e}^{- k\varepsilon \left(m+\lambda t\right)}}{k} $ , and the probability of “ultimate” transition to a new condition being $ \frac{1-{e}^{- k\varepsilon \left(m+\lambda t\right)}}{k} $ . We would accommodate more complex models of character state transitions such as those outlined earlier by varying ε in the same way that those models vary α. Note that we format equation (3) as we did because this is how it currently must be implemented to be used in programs such as RevBayes (Höhna et al. Reference Höhna, Landis, Heath, Boussau, Lartillot, Moore, Huelsenbeck and Ronquist2016; see Supplementary Code). The simplest depiction is:

(3a) $$ {Q}_P=\left(\begin{array}{cccc}-\varepsilon \left(m+\lambda t\right)& \varepsilon \left(\frac{m+\lambda t}{\left(k-1\right)}\right)& \dots & \varepsilon \left(\frac{m+\lambda t}{\left(k-1\right)}\right)\\ {}& & & \\ {}& & & \\ {}\varepsilon \left(\frac{m+\lambda t}{\left(k-1\right)}\right)& -\varepsilon \left(m+\lambda t\right)& \dots & \varepsilon \left(\frac{m+\lambda t}{\left(k-1\right)}\right)\\ {}& & & \\ {}\vdots & \vdots & \ddots & \vdots \\ {}& & & \\ {}\varepsilon \left(\frac{m+\lambda t}{\left(k-1\right)}\right)& \varepsilon \left(\frac{m+\lambda t}{\left(k-1\right)}\right)& \dots & -\varepsilon \left(m+\lambda t\right)\end{array}\right) $$

which reveals a major difference between the punctuated and continuous-time models that epitomizes the difference between punctuated equilibrium and “phyletic gradualism”: when t = 0 (i.e., a tree positing that a species is found immediately after a cladogenetic event or even a possible cladogenetic event), then we still expect some change. Indeed, if we somehow “knew” the true number of missing cladogenetic events along each branch, then time would cease to be a parameter here; however, as we cannot, we must use m + λt. Wagner and Marcot (Reference Wagner, Marcot, Alroy and Hunt2010) show that when t = 0, then the expected change converges on m + λt as we consider increasing numbers of possible events weighted by the probability of that many cladogenetic events given λt. Finally, note that it is not necessary to add a parameter to assess the probability that we have missed 1, 2, 3, etc., ancestors along each branch simply because this is already part of the prior probabilities of each branch assuming a fossilized birth–death or similar process (see Foote et al. Reference Foote, Hunter, Janis and Sepkoski1999; Heath et al. Reference Heath, Huelsenbeck and Stadler2014).

Cladogenesis Rates

Two issues immediately arise from tying expected character change to rates of cladogenesis (λ). One is whether allowing two parameters (λ and ε) instead of just one (α) to generate the likelihoods of possible evolutionary histories might unduly bias tests toward favoring punctuated equilibrium models. Another is how we go about including λ in phylogenetic analyses. Because the latter informs on the former, we will first review how Bayesian analyses incorporate diversification and sampling rates.

Bayesian phylogenetic analyses such as the fossilized birth–death (FBD) model (e.g., Stadler Reference Stadler2010; Heath et al. Reference Heath, Huelsenbeck and Stadler2014) assess the prior probability of branch durations on phylogenetic topologies based on rates of cladogenesis (λ), extinction (μ), and sampling (ψ) (Fig. 2). Sampling (ψ) alone affects the probability that we fail to sample either an earlier member of a lineage or some direct ancestor to one or more sampled taxa over time t (Wagner Reference Wagner1995a; Foote Reference Foote1996b). Cladogenesis (λ) affects the probability that 0…N cladogenetic events will happen over that time, and the combination of λ, μ, and ψ affect the probability that we would fail to sample any descendants of each cladogenetic event (Foote et al. Reference Foote, Hunter, Janis and Sepkoski1999; Stadler Reference Stadler2010; Bapst Reference Bapst2013; Didier et al. Reference Didier, Fau and Laurin2017). Thus, λ has a major effect on the posterior probability of a tree regardless of the model of character change (Fig. 2A), and using λ to co-determine expected character change does not introduce any new parameters into our model (Fig. 2B).

Figure 2. Graphical model for Bayesian analyses of evolutionary histories with varying rates of origination, extinction, and sampling over time and constant rates of change among characters and over time. In each graph, the likelihood component is on the left and is evaluated by character data, whereas the prior probability component is on the right and is evaluated by first appearance data. A, Continuous change, where the matrix of transition probabilities (Q) reflects the instantaneous rates of change (α) and duration of branches from a particular phylogeny (t). B, Punctuated change, where Q reflects the probability of change per branching event (ε) and t. Although allowing Q to reflect two parameters can increase the likelihood component of the posterior probability, this also demands that cladogenesis rate (λ) satisfies both first appearance and character data.

As described elsewhere (e.g., Wright Reference Wright2019; Wright et al. Reference Wright, Wagner and Wright2021, Reference Wright, Bapst, Barido-Sottani and Warnock2022), we can make numerous modifications to the basic models illustrated in Figure 2 that relax assumptions about homogeneity of rates. We would make identical modifications to both the continuous‑time and punctuated models in Figure 2 in order to allow for rate variation among characters (e.g., gamma or lognormal vs. invariant rates; e.g., Yang Reference Yang1994; Wagner Reference Wagner2012; Harrison and Larsson Reference Harrison and Larsson2015) or over the phylogeny (e.g., relaxed vs. strict clocks; e.g., Sanderson Reference Sanderson1997; Drummond et al. Reference Drummond, Ho, Phillips and Rambaut2006). Skyline models (e.g., Stadler Reference Stadler2011; Stadler et al. Reference Stadler, Kühnert, Bonhoeffer and Drummond2013; Warnock et al. Reference Warnock, Heath and Stadler2020) allow for the temporal variation in diversification (or sampling) rates of the sort that innumerable paleobiology studies document (e.g., Krug and Jablonski Reference Krug and Jablonski2012; Alroy Reference Alroy2015; Foote et al. Reference Foote, Cooper, Crampton and Sadler2018, Reference Foote, Edie and Jablonski2024; Lowery and Fraass Reference Lowery and Fraass2019; Foote Reference Foote2023). Under Skyline FBD models, an X Myr branch duration in an interval with high λ usually will have a lower probability than will an X Myr duration in an interval with low λ (see Wagner Reference Wagner2019). This is where the continuous-time and punctuated models might yield very different results, particularly if we are using “strict clock” models for both α and ε. Variation in λ over time will allow shifts in expected amounts of change per million years, even under constant ε. However, λ is simultaneously evaluated by the prior probabilities of branch durations: as λ increases, the probability of shorter branches increases, and as λ decreases, the probability of shorter branches increases (see, e.g., Wagner Reference Wagner2019). Thus, we cannot freely vary λ to maximize how well constant (i.e., strict clock) ε maximizes the probability of the character data. Instead, the punctuated model should result in high posterior probabilities relative to the continuous-time model only if increasing or decreasing origination rates to improve likelihood given character data also improves prior probabilities given stratigraphic data.

Data and Analyses

We examine 37 Ordovician genera from the brachiopod superfamily Strophomenoidea based on character data published by Congreve et al. (Reference Congreve, Krug and Patzkowsky2015; data available from Congreve et al. [Reference Congreve, Patzkowsky and Wagner2021b]) Strophomenoids appear during the Middle Ordovician and are one of the two superfamilies (along with the Plectambonitoidea) that contribute most to the diversification of “articulate” brachiopods during the Great Ordovician Biodiversification Event (GOBE) (Harper et al. Reference Harper, Rasmussen, Liljeroth, Blodgett, Candela, Jin, Percival, Rong, Villas and Zhan2013). Congreve et al. (Reference Congreve, Patzkowsky and Wagner2021a) show that the clade displays an early burst of morphological disparity that coincides with elevated rates of change per million years given phylogenies consistent with the cladogram published by Congreve et al. (Reference Congreve, Krug and Patzkowsky2015). However, those elevated rates themselves coincide with similar peaks in origination rates, leading Congreve et al. (Reference Congreve, Patzkowsky and Wagner2021a) to suggest that the early burst might reflect elevated frequencies of punctuated change rather than changes in the basic evolvability of strophomenoid shells. The former interpretation is more in keeping with predictions of “empty ecosytems” driving morphological change rather than “relaxed intrinsic constraints” (e.g., see Valentine Reference Valentine1969, Reference Valentine1980) and thus highlights the potential importance of being able to distinguish between continuous‑time and punctuated character change models when addressing other macroevolutionary issues.

We conduct analyses in RevBayes (Höhna et al. Reference Höhna, Landis, Heath, Boussau, Lartillot, Moore, Huelsenbeck and Ronquist2016). The following aspects of the punctuated and continuous‑time analyses are identical for both analyses. Both analyses assume lognormal rate variation among characters, as prior analyses indicate that this fits better than gamma or invariant distributions (Wagner Reference Wagner2012; Harrison and Larsson Reference Harrison and Larsson2015). Both use a 12-bin Skyline model, with chronostratigraphic partitions reflecting the Ordovician stage slices proposed by Bergström et al. (Reference Bergström, Chen, Gutiérrez–Marco and Dronov2009) and the Ordovician timescale given in Gradstein et al. (Reference Gradstein, Ogg, Schmitz and Ogg2020), although after lumping some of the relatively brief stage slices (see Congreve et al. Reference Congreve, Patzkowsky and Wagner2021a). We used sampling (ψ) and origination (λ) rates derived from published birth–death sampling analyses (Congreve et al. Reference Congreve, Patzkowsky and Wagner2021a) for the initial rates, but we then allowed these rates to vary over the analyses. These rates are given in the RevBayes scripts provided in the Supplementary Material. Following Wright et al. (Reference Wright, Wagner and Wright2021), we did not have extinction (μ) vary independently of origination; instead, we use a separate “turnover” parameter (i.e., μ/λ) that follows a lognormal distribution matching substage variation in relative extinction and origination rates found in empirical studies. Finally, the upper and lower bounds of the first appearance times for the analyzed taxa reflect the earliest and latest possible dates of the oldest species in each genus based on data in the Paleobiology Database (PBDB). However, instead of the dates given by the PBDB, we use an extrinsic database (see Congreve et al. Reference Congreve, Patzkowsky and Wagner2021a) that restricted possible first appearances to conodont or graptolite zones, with the ages of those zones again derived from Gradstein et al. (Reference Gradstein, Ogg, Schmitz and Ogg2020). We could not implement the complex compound beta distributions for possible first appearance dates used in Congreve et al. (Reference Congreve, Patzkowsky and Wagner2021a); therefore, the analyses assumed a uniform distribution of possible ages within those upper and lower bounds (see Barido-Sottani et al. Reference Barido-Sottani, Aguirre-Fernández, Hopkins, Stadler and Warnock2019). Barido-Sottani (Reference Barido-Sottani, Aguirre-Fernández, Hopkins, Stadler and Warnock2019; see also O’Reilly and Donoghue Reference O’Reilly and Donoghue2020) recommend using the entire stratigraphic range of taxa for the upper and lower bounds of a taxon’s stratigraphic range as the possible “true” first appearance age, noting that the failure to allow for this uncertainty biases the results. However, those authors did not consider “boundary crosser” taxa sensu Foote (Reference Foote2001) In this study, every analyzed taxon spans multiple conodont and graptolite zones, which means that only some of the finds might represent the oldest (or youngest) occurrences (see also Wright et al. Reference Wright, Wagner and Wright2021). For datasets like strophomenoid brachiopods, using the entire stratigraphic range biases results toward finding divergences after the oldest possible first appearance of a taxon.

For the purposes of this study, we use invariant rates (= strict clocks) for both continuous‑time (α) and punctuated (ε) rates of change. The sole difference in the two sets of analyses is that in the punctuated analysis, the rate of character change was set by ε, the “local” cladogenesis rate (λ _i) and the minimum branching events (m). For branches spanning 2+ stage slices, the expected change from cladogenesis in addition to m is given by:

(4) $$ \varepsilon \times \sum_{i=1}^n{\lambda}_i{t}_i $$

with ∑t_i equaling the branch duration. We accomplished this by: (1) dividing each branch into 100 kyr intervals; (2) referring to the origin time of each branch on each particular phylogeny; and (3) setting the rate of change to ελ _i where λ _i is the origination rate during that 100 kyr time slice (see Supplementary Material). We set m to 0.5 for all branches, because this was a genus-level analysis and we had no strong intuition of which branches the cladogenetic events should sit on. Even if a particular phylogeny has: (1) one genus as ancestral to another, and (2) both genera coexisting, it is possible that the branching event allowing 2+ lineages to coexist happened within the paraphyletic genus, while the “derived” genus represents greater amounts of continuous anagenetic change. This neglects the possibility that an ancestral genus might disappear before the first appearance of the derived genus. Unfortunately, current implementations of Markov chain Monte Carlo (MCMC) phylogenetic analyses do not include last appearance data and thus treat each taxon as a single occurrence in time (but see Warnock et al. Reference Warnock, Heath and Stadler2020). Fortunately, initial analyses (e.g., Congreve et al. Reference Congreve, Krug and Patzkowsky2015) do not suggest this happened within the clade, and most of the genera have wide ranges with apparent close relatives usually coexisting. As we will discuss later, using first and last appearance data that will allow phylogenetic algorithms to recognize the different scenarios illustrated in Figure 1 is an area that requires further development.

We set up each of our models in RevBayes 1.2.1 (Höhna et al. Reference Höhna, Landis, Heath, Boussau, Lartillot, Moore, Huelsenbeck and Ronquist2016). Additionally, we collected the following statistics from RevBayes: the Colless imbalance (Colless Reference Colless1995) of the trees, and tree length to be used for contrasting the results. We ran the models until convergence, as assessed in the software tracer (about 1 million generations per estimation). RevBayes returns both trees and parameter estimates in a tab‑separated format. We then processed results in the R programming environment using both our own scripts and existing programs (R Core Team 2018; Paradis and Schliep Reference Paradis and Schliep2019).

Ideally, we would contrast the results from the two models using tests such as stepping-stone analyses (Xie et al. Reference Xie, Lewis, Fan, Kuo and Che2011). Unfortunately, our implementation of the punctuated model results in very slow run times: whereas the continuous-time model requires approximately 12 hours on a fast laptop to execute 1 million MCMC generations, the punctuated model requires more than 2 months on the same device. Because methods such as stepping stone require multiple replications (usually 20), it was not possible to complete them using the punctuated model at the present time.