Inferindo variação da taxa de diversificação de filogenias com fósseis 

Jonathan S Mitchell, Rampal S Etienne, Daniel L Rabosky

Biologia Sistemática, Volume 68, Edição 1, Janeiro 2019, Páginas 1-18, https://doi.org/10.1093/sysbio/syy035

Publicado:

 

18 May 2018

 História do artigo

• PDF
•  Exibição dividida
•  Citar
• Ícone de permissões

   Permissões
•  
  Compartilhar

   

Abstrair

Filogenias calibradas por tempo de espécies vivas têm sido amplamente utilizadas para estudar o ritmo e o modo de diversificação de espécies. No entanto, está cada vez mais claro que inferências sobre a diversificação de espécies — taxas de extinção em particular — podem não ser confiáveis na ausência de dados paleontológicos. Introduzimos uma estrutura geral baseada no processo fossilizado de nascimento-morte para estudar a dinâmica de especiação - extinção em filogenias de espécies existentes e extintas. O modelo pressupõe que os filogenias podem ser modelados como uma mistura de regimes de taxa evolutiva distintos e que um processo hierárquico de Poisson rege o número de tais regimes de taxa em uma árvore. Implementamos o modelo no BAMM, uma estrutura computacional que usa a cadeia de salto reversível Markov Monte Carlo para simular uma distribuição posterior de regimes de taxa macroevolucionária condicionada aos tempos de ramificação e topologia de uma filogenia. A implementação, descrevemos pode ser aplicada a filogenias paleontológicas, filogenias neontológicas e a filogenias que incluem taxas extasiados e extintas. Avaliamos o desempenho do modelo em conjuntos de dados simulados em uma série de cenários de diversificação. Descobrimos que as taxas de especiação são inferidas de forma confiável na ausência de dados paleontológicos. No entanto, a inclusão de observações fósseis aumenta substancialmente a precisão das estimativas de taxa de extinção. Demonstramos que as inferências são relativamente robustas a pelo menos algumas violações de premissas de modelos, incluindo heterogeneidade nas taxas de preservação e a não especificação do número de ocorrências em conjuntos de dados paleontológicos.

BAMM, métodos comparativos, extinção, fósseis, macroevolução, especiação

Seção de problemas:

 Artigos regulares

Editor Associado: Richard Ree

A variação nas taxas de especiação e extinção contribuem para disparidades marcantes na riqueza de espécies entre diferentes clas de organismos. O estudo das taxas macroevolucionárias foi pioneiro no registro fóssil, com a aplicação de modelos fenomenológicos de nascimento-morte a padrões de diversidade ao longo do tempo (Raup et al. 1973; Sepkoski 1978; Raup 1985). Filogenias de espécies já existindo são amplamente utilizados para estimar taxas de especiação e extinção (Nee et al. 1994; Nee 2006; Alfaro et al. 2009; Morlon 2014). No entanto, a integração de dados paleontológicos e neontológicos provavelmente fornecerá insights mais precisos sobre a dinâmica macroevolucionária do que apenas os filogenias moleculares (Quental e Marshall 2010; Liow et al. 2010a; Stadler 2013; Rabosky 2016). De fato, um estudo de simulação recente argumentou que as taxas macroevolucionárias inferidas usando apenas uma topologia de árvores e datas de ocorrência fóssil são mais confiáveis do que as taxas inferidas com uma árvore e comprimentos de galhos apenas (Didier et al. 2017).

Inferências sobre padrões macroevolucionários frequentemente diferem entre dados filogenéticos paleontológicos e moleculares (Hunt e Slater 2016). Os dados fósseis suportam rotineiramente altos níveis de extinção e grandes flutuações na riqueza permanente (por exemplo, Alroy 1999, 2008, 2010; Liow et al. 2010a; Quental e Marshall 2010; Folhas et al. 2016). No entanto, análises somente extasiais raramente encontram suporte para altas taxas de extinção (Nee 2006; Purvis 2008, mas veja Etienne et al. 2012), e inúmeros estudos abordaram possíveis razões para essa discrepância (Rabosky 2009, 2010, 2016; Morlon et al. 2011; Etienne e Rosindell 2012; Etienne et al. 2012; Rosenblum et al. 2012; Beaulieu e O'Meara 2015). A incorporação direta de dados fósseis na estimativa da taxa filogenética deve facilitar testes mais robustos do papel da extinção na dinâmica macroevolucionária (Sakamoto et al. 2016). Uma abordagem explicitamente filogenética para taxas macroevolucionárias no registro fóssil fornecerá insights sobre a variação da taxa entre as linhagens e também facilitará o estudo da dinâmica de diversificação em clados que são muito mal preservados para abordagens baseadas em subsampling (por exemplo, Alroy 1999).

Neste artigo, descrevemos uma estrutura geral para modelar dinâmicas complexas de especiação - extinção em árvores filogenéticas que incluem dados paleontológicos e neontológicos. O modelo é baseado no processo fossilizado de nascimento e morte (Stadler 2010; Didier et al. 2012, 2017; Heath et al. 2014) e inclui parâmetros para especiação, extinção e taxas de fossilização. O processo de nascimento-morte para árvores filogenéticas reconstruídas (Nee et al. 1994) é um caso especial do processo fossilizado de nascimento-morte onde a taxa de preservação é igual a zero (Stadler 2010). Como tal, este modelo pode ser aplicado a filogenias paleontológicas (apenas taxa extinta), filogenias neontológicas de espécies vivas (somente taxa de subtaranidade) e filogenias que incluem qualquer combinação de espécies vivas e extintas.

A nova característica de nossa implementação, em relação ao trabalho anterior sobre taxas de diversificação de filogenias e fósseis (Stadler 2010; Didier et al. 2012; Heath et al. 2014), é a capacidade de capturar padrões complexos de variação entre linhagem nas taxas de especiação e extinção. O modelo pressupõe que os filogenias são moldados por um conjunto contavelmente finito de regimes de taxa evolutiva distintos, e que o número de tais regimes é regido por um processo hierárquico de Poisson. Este modelo é uma extensão do modelo de diversificação multiprocesso atualmente implementado no pacote de software BAMM para análise de taxas macroevolucionárias (Rabosky 2014; Rabosky et al. 2017). Notamos que Sakamoto et al. (2016) desenvolveram um quadro flexível baseado em regressão para estimar a heterogeneidade nas taxas de diversificação de filogenias fósseis, mas sua abordagem não se baseia em um modelo formal baseado em processos de diversificação e preservação.

Primeiro descrevemos a formulação analítica do modelo e afirmamos suas suposições. Em seguida, descrevemos nossa implementação dentro da estrutura de software BAMM e discutimos sua relação com as versões BAMM anteriores (Rabosky et al. 2017). Para avaliar o desempenho do modelo, simulamos filogenias sob uma série de cenários de diversificação e, em seguida, submetemos cada filogenia simulada à fossilização estocástica para gerar uma árvore semelhante ao que os pesquisadores normalmente usariam (uma árvore filogenética "observada" com uma história incompletamente fossilizada). Analisamos cada filogenia com BAMM e comparamos as taxas inferidas pelo BAMM e as configurações de mudança aos valores conhecidos (simulados) para avaliar a precisão das estimativas dos parâmetros. Descobrimos que a inclusão de dados paleontológicos até mesmo esparsos em filogenias dominadas por linhagens já controladas melhora a precisão das estimativas de especiação e taxa de extinção. Essa melhora é especialmente acentuada na estimativa das taxas de extinção. Embora nosso modelo faça várias suposições simplificadoras sobre a natureza da fossilização, demonstramos que as inferências de taxa são razoavelmente robustas à heterogeneidade temporal nas taxas de preservação.

Métodos

Dados

Os dados consistem em uma árvore filogenética calibrada pelo tempo (que pode incluir linhagens extintas e extasiadas) e uma contagem do número de ocorrências fósseis estratigraficamente distintas associadas à árvore. Cada linhagem extinta na filogenia deve ser representada por um ramo que termina no momento da última ocorrência para o taxon. Observe que os verdadeiros eventos de extinção não são observados diretamente, e espera-se que o último tempo de ocorrência preceda o verdadeiro tempo de extinção de um taxon. O número de ocorrências fósseis é melhor pensado como uma estimativa do número de fósseis que poderiam ser, em princípio, atribuídos a linhagens presentes na árvore. Ocorrências fósseis de natureza indeterminada e potencialmente associadas a linhagens não presentes na árvore (por exemplo, "Theropoda, indeterminada") não são acomodadas neste quadro de modelagem e são excluídas da contagem de ocorrências; os cálculos de probabilidade só se aplicam a fósseis que podem ser atribuídos explicitamente a galhos na árvore. Dito de outra forma, as ocorrências fósseis só devem ser contadas se o imposto de nível de espécie a que são atribuídos for incluído na árvore. No entanto, como a implementação atual pressupõe homogeneidade de preservação através do tempo e entre as linhagens, os cálculos não utilizam de fato contexto estratigráfico e filogenético para observações fósseis individuais. Como explicado abaixo, essa suposição de preservação homogênea do tempo significa que a probabilidade da filogenia com dados fósseis é invariante em relação à localização dos fósseis na árvore.

Likelihood of a Phylogeny with Fossils and Rate-Shifts

The fundamental operation in the analysis of diversification rate shifts on phylogenetic trees involves computing the likelihood of the phylogeny under a fixed configuration of lineage types. Each “type” of lineage is a characterized by a distinct and potentially time-varying rate of speciation (⁠

λ)

$\lambda )$ and extinction (⁠μ)

$\mu )$⁠. Two lineages i

$i$ and 

j

$j$ are said to be of the same type if λi

(t)=

${\lambda }_i(t)=$λ

j

(t)

${\lambda }_j(t)$ and μi

(t)=

μ

j

(t)

${\mu }_i(t)={\mu }_j(t)$ everywhere in time. Thus, for a given point in time 

t

$t$⁠, lineages of the same type will have precisely the same distribution of progeny lineages at some time 

t+Δt

$t+\mathrm{\Delta }t$ in the future. For the constant-rate birth–death process, all lineages are assumed to be of the same type. For the Binary State Speciation and Extinction (BiSSE) process (Maddison et al. 2007), phylogenies are assumed to comprise a mixture of two types of lineages, and the character states at the tips of the tree can be viewed as labels for lineage types. The transition points between types are unknown in the BiSSE process, and the likelihood of a given set of labeled types can be computed by integrating over all possible transitions that could have given rise to the observed data.

In contrast to BiSSE, the calculations implemented in BAMM and other rate-shift models (Alfaro et al. 2009; Morlon et al. 2011; Etienne and Haegeman 2012) assume that we have mapped lineage types across all branches of the phylogeny and not merely to the tips. We define a “mapping” of rate parameters as an association between each point on a phylogeny (e.g., segment of a branch) and a particular lineage type. Additional steps may be used to optimize the number, location, and parameters of the mapped types, including maximum likelihood (Alfaro et al. 2009; Morlon et al. 2011; Etienne and Haegeman 2012).

Likelihood calculations for the fossilized birth–death process were described by Stadler (2010), and our general approach extends these calculations to a phylogeny with multiple types of lineages. The fossilized process differs from the reconstructed evolutionary process (Nee et al. 1994) because we must also account for the lineage fossilization rate 

ψ

$\psi$⁠. The new calculations described below are implemented in BAMM v2.6.

The derivation of the likelihood for rate-shift models follows the logic described in Maddison et al.’s (2007) description of the BiSSE model. In a given infinitesimal interval of time 

Δt

$\mathrm{\Delta }t$⁠, the state-space for the process includes five events that can potentially change the probability of the data: a lineage may undergo a rate shift, become preserved as fossil, become extinct, speciate, or it may undergo none of these events. To compute the likelihood of the tree and a set of parameters (e.g., a mapping of rate shifts), we assume that no other rate shifts have occurred, thus dropping one of the possible events (the occurrence of an unmapped rate shift) from the state space for the process. As such, the likelihood is conditioned on the set of shifts that have been placed on the tree. This assumption of rate-shift models has been criticized in the recent literature (Moore et al. 2016). However, it is unclear how we might accommodate unobserved rate shifts, even in principle, as there is very little information in the data or elsewhere that can yield information about the probability distribution of unobserved rate parameters. Most importantly, we demonstrate that the simplifying assumptions that underlie the BAMM calculations yield robust inference on evolutionary rates, even when simulated phylogenies used for validation contain many such unobserved rate shifts.

We assume that preservation events are mutually exclusive with respect to speciation and extinction events. Hence, a lineage cannot leave a recoverable fossil observation during the precise moment of lineage splitting or lineage extinction. This interpretation of the preservation rate ensures that the model we have implemented is mathematically identical to the model described by Stadler (2010), for the special case where speciation and extinction rates do not vary among lineages. This interpretation of preservation is also consistent with paleontological practice: fossil observations are always assigned to single lineages. Even in the relatively high-resolution marine microfossil record, paleontologists do not assume that the data capture the precise instant of lineage splitting (Benton and Pearson 2001). In this article, we assume that lineage preservation is a time-homogeneous process, but it is straightforward to relax this assumption.

We measure time 

t

$t$ backwards from some arbitrary starting time, 

t

0

$t_0$⁠. The start time, 

t

0

$t_0$⁠, need not correspond to the present. Let D(t)

$D(t)$ be the probability density that a lineage that is extant some time 

t

$t$ before 

t

0

$t_0$ gives rise to the observed data at the present time, and let E(t)

$E(t)$ be the corresponding probability that a lineage that is extant at time 

t

$t$ goes extinct, along with its descendants, before time 

t

0

$t_0$⁠. Following Stadler (2010), the master equations governing the likelihood of the data are

dD

dt

=−(λ+μ+ψ)D(t)+2λD(t)E(t)

$$\frac{dD}{dt}=-\left(\lambda +\mu +\psi \right)D\left(t\right)+2\lambda D\left(t\right)E\left(t\right)$$

(1)

and

dE

dt

=μ−(λ+μ+ψ)E(t)+λE

(t)2

$$\frac{dE}{dt}=\mu -\left(\lambda +\mu +\psi \right)E\left(t\right)+\lambda E{\left(t\right)}^2$$

(2)

E(t)

$E(t)$ describes the probability that an independent lineage originating at time 

t

$t$ goes extinct, along with all of its descendants, before some future time 

t

0

$t_0$⁠. One can also view this probability as the probability of “non-observation” of a lineage and all of its descendants; a lineage and its descendants might not be observed because they have gone extinct, and/or because we have failed to sample them. To compute the likelihood of the phylogeny with mapped rate shifts, we require solutions to these equations for arbitrary 

t

$t$ and initial values 

t

0

$t_0$⁠, D0

$D_0$⁠, and 

E

0

$E_0$⁠. For E(t)

$E(t)$⁠, the solution is (Stadler 2010)

E(t)=

λ+μ+ψ+

c1e

−

c1

Δt

(1−

c2

)−(1+

c2

)

e

−

c1

Δt

(1−

c2

)+(1+

c2

)

2λ

$$E\left(t\right)=\frac{\lambda +\mu +\psi +c_1\frac{e^{-c_1\mathrm{\Delta }t}\left(1-c_2\right)-\left(1+c_2\right)}{e^{-c_1\mathrm{\Delta }t}\left(1-c_2\right)+\left(1+c_2\right)}}{2\lambda }$$

(3)

where

c1

=

∣∣∣(λ−μ−ψ)2

+4λψ

‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾

√∣∣∣

$$c_1=\left|\sqrt{{\left(\lambda -\mu -\psi \right)}^2+4\lambda \psi }\right|$$

(4)

and

c2

=−

λ−μ−ψ−2λ(1−

E

0

)

c1

.

$$c_2=-\frac{\lambda -\mu -\psi -2\lambda \left(1-E_0\right)}{c_1}.$$

(5)

For 

D(t

$D(t$⁠, 

t

0

)

$t_0)$⁠, we have

D(t)=

4

D0

2(1−

c22

)+

e

−

c1

Δt

(1−

c2

)

2

+

ec1

Δt

(1+

c2

)

2

$$D\left(t\right)=\frac{4D_0}{2\left(1-c_2^2\right)+e^{-c_1\mathrm{\Delta }t}{\left(1-c_2\right)}^2+e^{c_1\mathrm{\Delta }t}{\left(1+c_2\right)}^2}$$

(6)

which follows immediately from Stadler (2010) Theorem 3.1, under the substitution 

D0

=ρ

$D_0=\rho$⁠. In the Supplementary Material, we provide a derivation of D(t)

$D(t)$ for arbitrary D0

$D_0$⁠.

Calculations are performed by solving 

D(t)

$D(t)$ and E(t)

$E(t)$ recursively on the phylogeny from the tips to the root, combining D(t)

$D(t)$ probabilities at internal nodes (discussed below) and continuing down the tree in a rootwards direction. For lineages that are extant, we have initial conditions at time 

t

0

=

$t_0=$ 0 of D0

=ρ

$D_0=\rho$ and 

E

0

=

$E_0=$ 1−ρ

$-\rho$⁠, where ρ

$\rho$ is the probability that an extant lineage has been sampled (e.g., the sampling fraction; FitzJohn et al. 2009). In the absence of perfect preservation, the last fossil occurrence of a given taxon is not the true extinction time. For a lineage with a last occurrence at time 

t

i

$t_i$⁠, we generally know that the lineage and all of its descendants went extinct sometime between 

t

i

$t_i$ and 

t

0

=

$t_0=$ 0. Hence, for terminal lineages with non-extant tips, we initialize both 

E

0

$E_0$ and D0

$D_0$ with 

E(

t

i

∣

t

0

=

$E(t_i\mid t_0=$ 0), or the probability that a single lineage that is extant at time 

t

i

$t_i$ will have zero observed descendants prior to and including the present day (⁠t=

$t=$ 0). This initialization for D0

$D_0$ accounts for the probability that a non-extant tip goes extinct, or is unsampled, on the time interval (⁠

t

i

$t_i$⁠, 0). At internal nodes, we combine the densities from the right DR

(t)

$D_R(t)$ and left D

L

(t)

$D_L(t)$ descendant lineages, such that D(t)

$D(t)$ immediately rootwards of the node is computed as 

D(t)=λ(t)

DR

(t)

D

L

(t)

$D(t)=\lambda (t)D_R(t)D_L(t)$⁠. The full likelihood must also account for all 

Z

$Z$ observed fossil occurrences, which is accomplished by multiplying the probability density of the tree by ΨZ

${\mathrm{\Psi }}^{\mathrm{Z}}$⁠.

There is one difference between the implementation of the model described here (implemented in BAMM v2.6) and the default calculation scheme implemented in the most recent previous version, BAMM v2.5. As discussed in Rabosky et al. (2017), there are several algorithms for computing the likelihood of a specified mapping of rate regimes on a phylogenetic tree. These algorithms, which we have termed “recompute” and “pass-up” (Rabosky et al. 2017), differ in how the analytical 

E(t)

$E(t)$ equation is initialized for calculations on individual branch segments. The precise sequence of calculations associated with these approaches is described in the Supporting Information (Section S2.4 available on Dryad at https://doi.org/10.5061/dryad.50m70) from Rabosky et al. (2017). Herein, we expand BAMM to include the “recompute” algorithm for computing the likelihood, which was not implemented as the default calculation scheme in BAMM v2.5. The process under consideration in the present article can yield phylogenies that have gone extinct in their entirety (e.g., fossil taxa only), and we can thus avoid conditioning on survival of the process as well as unresolved theoretical problems associated with this conditioning. We thus perform all calculations using both the “recompute” and “pass-up” algorithms; we present “recompute” results in the main text, but the corresponding results using “pass-up” were virtually identical and are presented in the Supplementary Information available on Dryad.

Como observado acima, nosso uso dessas equações para modelar dinâmicas de diversificação condie implicitamente a probabilidade em um conjunto finito e observado de tipos de linhagem: a probabilidade não explica novos tipos de linhagens que podem ter evoluído em linhagens sem descendentes fossilizados ou existentes (Moore et al. 2016). Essa suposição se aplica a todos os métodos que pretendem calcular a probabilidade de uma filogenia sob um mapeamento especificado de mudanças de taxa (Alfaro et al. 2009, Morlon et al. 2011, FitzJohn 2010, Etienne e Haegeman 2012). De fato, essa suposição se aplica indiretamente a todos os métodos de diversificação: se não houver (ou muito poucas) mudanças não observadas, a probabilidade calculada pelo BAMM é aproximadamente correta. Se a "realidade" contiver muitas dessas mudanças não observadas, então a probabilidade (e as estimativas de taxa correspondentes) do BAMM serão tendenciosas, porque as mudanças de taxa terão essencialmente servido como uma fonte não acomodada de extinção da linhagem. Mas se mudanças de taxa não observadas são suficientemente comuns (em dados reais) quanto ao viés BAMM, então elas são problemáticas para todos os outros métodos e modelos atuais também. Como apontado por Rabosky et al. (2017), o problema das mudanças de taxa não observadas não desaparece simplesmente porque se assume um modelo teoricamente completo que ignora sua existência ou por amostragem de seus parâmetros de uma distribuição anterior (assumida).

Implementação no BAMM

O modelo descrito acima é implementado como uma extensão ao software BAMM (v2.6). O BAMM v2.6 expande o modelo básico de diversificação implementado em versões anteriores do BAMM, incorporando parâmetros que especificam o tempo total de observação da árvore ("observaçãoTime", 

T

obs

${}_{\mathrm{Tobs}}$⁠: tempo da raiz quando a árvore é observada; este parâmetro é igual à idade raiz das árvores com linhagens existentes), o número total de ocorrências fósseis observadas dentro do clado analisado ("numerações ocorrências"), e parâmetros que controlam a distribuição prévia e atualização da frequência para a taxa de preservação Ψ

$\mathrm{\Psi }$ ("updatePreservationRate", "preservationRatePrior" e "preservationRateInit").

T

obs

${}_{\mathrm{Tobs}}$ pode ter um grande impacto no ajuste geral do modelo, e garante uma explicação adicional. O modelo pressupõe que as linhagens de ponta não-excidente podem ser extintas em qualquer lugar no intervalo de tempo que se estende desde sua última ocorrência até 

T

obs

${}_{\mathrm{Tobs}}$⁠. Considere o cenário em que estamos modelando dinâmicas de diversificação em um clado extinto de dinossauros que se originou 200 milhões de anos antes do presente (Ma) e onde a última ocorrência desse clado é de 90 Ma. Se definirmos 

T

obs

${}_{\mathrm{Tobs}}$ para igualar os dias atuais (0 Ma), o modelo permitirá a possibilidade de que o último dinossauro do clade foi extinto em algum momento entre 90 e 0 Ma. Na prática, sugerimos a escolha da observaçãoTime como o mais cedo possível que acredita-se que todos os impostos de última ocorrência tenham sido extintos. No exemplo do dinossauro, um cenário razoável para 

T

obs

${}_{\mathrm{Tobs}}$ seria o limite de KPg em 66 Ma, quando todas as linhagens de dinossauros não-aviários haviam se extinto.

O número de ocorrências (parâmetro "numerações" no BAMM) representa o número total de ocorrências fósseis que podem ser atribuídas a linhagens observadas no clade em consideração e não é simplesmente o número de pontas extintas na filogenia. É necessariamente o caso de que um filogenia com 

k

$k$ dicas extintas tem, pelo menos, 

k

$k$ocorrências: uma para a última ocorrência de cada linhagem extinta; no entanto, pode haver múltiplas ocorrências para cada linhagem extinta e exciente. Existem várias interpretações possíveis de uma ocorrência fóssil no quadro BAMM; abordamos este assunto longamente na seção Discussão. As ocorrências associadas à taxa fóssil que não estão presentes na árvore não devem ser utilizadas, nem material taxonomicamente indeterminado (por exemplo, material fóssil fragmentário que não pode ser identificado taxonomicamente abaixo do nível do gênero). Em geral, consideramos uma "ocorrência" para representar observações estratigraficamente únicas em nível de espécie. Assim, uma localidade com 1000 indivíduos de um único taxon representando um único evento deposição (por exemplo, um evento de enterro em massa) constituiria uma única ocorrência. Na prática, sugerimos analisar um conjunto de dados com vários valores diferentes para a contagem de ocorrências para garantir que os resultados sejam robustos (ver seção Discussão).

Descrição da Simulação de Árvores

Desenvolvemos software para simular filogenias sob um processo poisson que inclui mudanças para novos regimes de taxa macroevolucionária (https://github.com/macroevolution/simtree). Cada simulação foi iniciada com duas linhagens ao mesmo tempo 

t=

$t=$ 0 com λ

$\lambda$ e μ

$\mu$ extraídos de distribuições exponenciais com médias de 0,2 e 0,18, respectivamente. Assumimos que mudanças para novos regimes de taxa ocorreram com taxa Λ=

$\mathrm{\Lambda }=$ 0.02. Waiting times between events (speciation, extinction, rate shift) were drawn from an exponential distribution with a rate equal to the sum of the rate parameters (⁠λ+μ+Λ)

$\lambda +\mu +\mathrm{\Lambda })$⁠. For rate shift events, new values of λ

$\lambda$ and μ

$\mu$ were drawn from the same exponential distributions as the root regime. This procedure simulates a phylogeny with stochastic rate shifts; within a given rate regime, λ

$\lambda$ and μ

$\mu$ were constrained to be constant in time. Simulations were run for 100 time units, and trees with more than 5000 or fewer than 50 lineages at t=

$t=$ 100 were discarded, resulting in a total of 200 trees for analysis. In addition, each tree was constrained to have at least one rate shift in the true, generating history. This simulation procedure can generate phylogenies with unobserved rate shifts when shifts occur in unsampled extinct clades (see below). Conversely, the true root regime may be unobserved and all observed tips may fall within one of the derived regimes. As discussed in the section below on fossilization, any rate shift that leads to an extinct clade with no fossil observations will be unobserved, thus enabling us to assess whether Moore et al.’s (2016) concerns about unobserved rate shifts have consequences for inference in practice.

Description of Fossilization Procedure

We created pseudo-fossilized trees by stochastically placing fossils along the branches of each simulated tree (Fig. 1) using per-lineage preservation rates (⁠

ψ)

$\psi )$ of 0, 0.01, 0.1, or 1.0. Speciation rates in our simulations ranged from 0.003 to 1.7 lineages per unit time, with 99% of the regimes having a speciation rate above 0.015; 80% exceeded 0.1. As such, speciation events were generally more frequent than fossilization events, except under the high preservation (⁠ψ=

$\psi =$ 1.0) simulation scenario. All extant tips, and tips with at least one fossil occurrence, were retained while tips (and clades) lacking fossil occurrences were dropped. In empirical data sets, the true extinction time of an extinct lineage is unknown, and the time of the last (most recent) fossil occurrence is used. For comparability with empirical data from the fossil record, extinct tips with fossils were truncated to the time of the most recent fossil occurrence (e.g., the branch ends at the last occurrence, and any subsequent history for the lineage is necessarily unobserved). We refer to these pseudo-fossilized trees as “degraded,” because the stochastic fossilization history necessarily discards information about the true lineage history.

Figure 1.

Open in new tabDownload slide

Diagram showing a complete evolutionary tree with unobserved (gray) and observed (black) evolutionary history. Fossils are indicated by open circles; labeled circles (a, b) denote rate shifts on individual branches. Tree contains three extant tips, three observed extinct tips, and eight unobserved extinct tips. Extinct taxa are necessarily placed at the time of their last fossil occurrence. Rate shift (a) is unobserved in phylogenies that include fossil data, as there are no extant descendants or fossil occurrences on the subtree descended from the shift. The simulation protocol used in this study involved generating phylogenies for the complete process, then simulating fossil occurrences on the tree and removing unobserved (gray) evolutionary history.

As pointed out by Moore et al. (2016), the model implemented in BAMM assumes that rate shifts do not occur on unobserved branches. However, our simulation and fossilization procedure allows rate shifts to occur on extinct clades with no fossil observations. Hence, our fossilization procedure frequently leads to unobserved rate regimes and enables us to test the practical significance of this core BAMM assumption.

BAMM Analyses

We approximated the posterior distribution of rate shift configurations for each degraded tree by running a Markov chain Monte Carlo (MCMC) simulation in BAMM v2.6 for 50 million generations, with the first 5 million discarded as burnin. Each tree was analyzed assuming a single, time-constant rate of preservation. The model prior (mean of the geometric prior distribution on the number of shifts) was set to 200 for all runs (after Mitchell and Rabosky 2016), while the rates of the exponential priors on 

λ

$\lambda$ and μ

$\mu$ were scaled using the setBAMMpriors() function in BAMMtools (see Supplementary Material available on Dryad for details). Mitchell and Rabosky (2016) found that inferences on the number of shifts were generally robust to the choice of model prior, but that use of a liberal model prior facilitated convergence of the MCMC simulation. We performed each analysis using both “recompute” and “pass-up” algorithms for computing the likelihood, as described by Rabosky et al. (2017). All computer code, data files, and result files are available through the Dryad digital data repository (https://doi.org/10.5061/dryad.50m70).

Performance Assessment

To assess the performance of our model and BAMM implementation, we quantified four fundamental attributes of interest: 1) the estimated preservation rate, 2) the speciation and extinction estimates inferred for each (true) rate regime, 3) the inferred number of macroevolutionary regimes on the tree, and 4) the location of inferred rate shifts. We compared the posterior distribution of the estimated preservation rate 

ψ

$\psi$ for each degraded tree to the generating value. For diversification rates we found the mean inferred speciation and extinction rates for every true regime with at least two tips (both extinct and extant), and compared these mean values to the true values of speciation (⁠λ)

$\lambda )$⁠, extinction (⁠μ)

$\mu )$ and relative extinction (⁠μ

$\mu$/λ)

$\lambda )$⁠. We predicted that BAMM would estimate rates more accurately for regimes with larger numbers of taxa; we thus assessed BAMM’s accuracy with explicit reference to the number of tips in each regime. This allowed us to assess how well BAMM estimates rates for regimes that differ in size.

We determined the best-supported number of shifts for each tree by using stepwise model selection with Bayes factors (Mitchell and Rabosky 2016). Using the posterior distribution of shifts for each tree, we compared the lowest complexity model (e.g., smallest number of shifts) to the next-most complex model. If the more complex model was supported by Bayes factor evidence of at least 20 relative to the simpler one, it was provisionally accepted and the process repeated with models of increasing complexity until the Bayes factor threshold of 20 was no longer exceeded. We then compared the best-supported number of shifts on the tree to the true number of preserved shifts (see Fig. 1).

Finally, we looked at how accurately BAMM inferred the locations of shifts along the tree. Location accuracy is difficult to quantify, as many simple metrics have undesirable properties. For example, the distance along the branches between a true and inferred shift is difficult to interpret unless the number of inferred shifts is exactly equal to the number of true shifts. Branch lengths further confound location accuracy as shifts along shorter branches are expected to be more “smeared” across adjacent branches relative to shifts on longer branches (see extensive discussion of this topic in the supplement to Mitchell and Rabosky 2016). We thus assessed location accuracy using several complementary metrics.

The first metric we used for location accuracy was the marginal posterior probability of at least one rate shift on all branches that included a (true) rate shift. We assessed these marginal probabilities with respect to both the magnitude of the shift (e.g., size of the rate increase) as well as the number of tips descending from that shift, with the expectation that shifts with many descendent tips and/or with a large change in rates would have higher marginal posterior probabilities. As pointed out by Rabosky (2014), focusing on these marginal probabilities can be misleading if the evidence for a rate shift is “smeared” across a set of adjacent branches.

As an alternative metric of location accuracy, we used the method described by Mitchell and Rabosky (2016) to assess how well BAMM partitions tips into their correct rate classes. This approach can be motivated by noting that any placement of 

k

$k$ rate shifts on a phylogeny partitions a tree into at most k+

$k+$ 1 subsets of tips with distinct rate regimes (⁠k+

$k+$ 1, not simply 

k

$k$⁠, due to the presence of a root regime that is distinct from the rate shifts). For a phylogeny with 

N

$N$ taxa and a known (true) mapping of diversification regimes, there is an 

N

$N$ x 

N

$N$ matrix where each element 

C

i,k

$C_{i,k}$ describes whether the i

$i$’th and 

k

$k$’th taxa are assigned to the same or different rate regimes (the set of tips partitioned into a particular rate regime was termed a “macroevolutionary cohort” in Rabosky 2014). Elements of the matrix 

C

i,k

$C_{i,k}$ take a value of 1 (if taxa i

$i$ and 

k

$k$ are in the same rate regime), or 0 (if i

$i$ and 

k

$k$ are in different rate regimes). Following BAMM analysis, we obtain a matrix D

$D$ where D

i,k

$D_{i,k}$ represents the posterior probability that taxa i

$i$ and 

k

$k$ are assigned to the same rate regime. To assess shift location accuracy, we computed one minus the average of the absolute value of the difference between the true pairwise partitioning matrix 

C

$C$ and the BAMM-estimated matrix D

$D$⁠, or

1−

2N(N−1)

∑k=2

N

∑i=1k−1∣∣

C

i,k

−

D

i,k

∣∣

$$1-\frac{2}{N\left(N-1\right)}\sum _{k=2}^N\sum _{i=1}^{k-1}\left|C_{i,k}-D_{i,k}\right|$$

(7)

An overall value of 1.0 can only be obtained if all pairs of taxa are correctly partitioned in the correct configuration of rate regimes in 100% of samples from the posterior. This metric of location accuracy will penalize shift configurations that include too many shifts (by separating taxa that should be in the same regime) and that include too few shifts (by uniting taxa that should be in separate regimes). However, the topology of the tree may make certain partitions more likely by chance alone. To quantify this baseline accuracy, we computed the partitioning accuracy that would be expected by chance if shifts were distributed randomly throughout the tree under a uniform prior with respect to the total tree length. Under a uniform prior, the probability of a shift on a given branch 

bi

$b_i$ is simply bi

$b_i$/

T

$T$⁠, where bi

$b_i$ is the length of the i

$i$’th branch and 

T

$T$ is the tree length. We conditioned the randomization on the simulated shift count, such that each randomized sample contained exactly the same number of shifts as the corresponding sample from the posterior. After randomizing the placements of shifts from all samples in the posterior, we computed the accuracy of random partitions for each sampled generation using the same pairwise distance matrix approach described above.

Sensitivity of the Method to Model Violations

We tested sensitivity of the framework described here to four qualitatively distinct violations of the assumptions of the underlying inference model. These violations included: 1) violation of time-homogeneous preservation potential; 2) misspecification of the true number of fossil occurrences in the data; 3) presence of a mass extinction in the data; and 4) diversity-dependent diversification. In real data sets, preservation rates can vary substantially even over relatively short timescales (Foote and Raup 1996; Holland 2016). We violated the assumption of time-constant preservation by applying a “stage-specific” preservation pattern to the simulated phylogenies. We use “stage” in this context solely to denote an arbitrary temporal span with a potentially distinct preservation rate. We divided each simulated tree into ten equally-sized temporal bins (each ten time units long, mirroring the 10 myr bins used in some large-scale paleodiversity studies), where each bin was assigned an independent preservation rate drawn from a uniform (0.01, 0.9) distribution. We then simulated fossilization histories using these preservation rates and pruned unobserved lineage histories as described above to generate a set of degraded trees. We analyzed these trees using BAMM v2.6, which assumes time-homogeneous preservation, and we tested the accuracy of inference using the same metrics described above. Given that the simulated trees involve rate shifts, this simulation approach allows clades to radiate in periods of higher or lower preservation, and to persist into subsequent time bins with distinct and decoupled preservation parameters.

Our implementation assumes that users can specify a meaningful number of fossil occurrences for the clade. Fossil occurrences are not placed at specific points on the tree, but rather inform the model about the total number of fossil occurrences associated with the focal clade. The current implementation thus assumes fossil occurrences (with the exception of last occurrences) are distributed across the tree under a homogeneous Poisson process, and estimates a tree-wide rate accordingly. However, for real data sets, the “true” number of occurrences is likely impossible to know with precision (even in principle; see Discussion section). To assess how misspecification of the number of fossil occurrences may bias inference, we simulated 100 constant rate phylogenies, with higher rates than the variable-rate tree simulations (mean 

λ=

$\lambda =$ 0.9, mean μ=

$\mu =$ 0.45) and imposed a fossilization history with 

ψ

$\psi$ of 0.1. Given the stochastic nature of the fossilization simulation, the fossilized trees varied in the number of extinct lineages observed. Note that because an extinct tip is only observed if at least one fossil occurs on that branch, each fossilized tree with (observed) extinct lineages has a minimum number of fossil occurrences equal to the number of (observed) extinct tips (FMIN

)

${}_{\mathrm{M}\mathrm{I}\mathrm{N}})$⁠. We then calculated the difference between this per-tree minimum and the true number (FTRUE

${}_{\mathrm{T}\mathrm{R}\mathrm{U}\mathrm{E}}$−

$-$ FMIN

=

${}_{\mathrm{M}\mathrm{I}\mathrm{N}}=$ FGAP

)

${}_{\mathrm{G}\mathrm{A}\mathrm{P}})$⁠, and analyzed the tree using BAMM v2.6 with the observed number of occurrences set to: 1) FMIN

${}_{\mathrm{M}\mathrm{I}\mathrm{N}}$⁠, 2) FMIN

+

${}_{\mathrm{M}\mathrm{I}\mathrm{N}}+$ 0.25 * FGAP

${}_{\mathrm{G}\mathrm{A}\mathrm{P}}$⁠, 3) FMIN

+

${}_{\mathrm{M}\mathrm{I}\mathrm{N}}+$ 0.5 * FGAP

${}_{\mathrm{G}\mathrm{A}\mathrm{P}}$⁠, 4) FTRUE

${}_{\mathrm{T}\mathrm{R}\mathrm{U}\mathrm{E}}$⁠, 5) FTRUE

+

${}_{\mathrm{T}\mathrm{R}\mathrm{U}\mathrm{E}}+$ FGAP

${}_{\mathrm{G}\mathrm{A}\mathrm{P}}$⁠, and 6) 10 * FTRUE

${}_{\mathrm{T}\mathrm{R}\mathrm{U}\mathrm{E}}$⁠. Overestimation of the true number of fossils may seem unlikely, but this could arise in practice due to taxonomic error (e.g., incorrect assignment of fossil observations to a specific taxon) or to the inclusion of indeterminate fossil material in the analysis. Scaling the degree of misspecification to the true and minimum number of fossils created roughly comparable degrees of misspecification across trees with large differences in the number of observed extinct lineages.

We estimated the bias for each of these misspecified fossil occurrence values by dividing the inferred speciation and extinction rates across the entire tree (as inferred using the getCladeRates function from BAMMtools) by the true values used to simulate the tree. Our discussion below pertains to bias in rates resulting from misspecification of the number of occurrences, not bias due a particular magnitude of preservation (e.g., “low” or “high” preservation). We predicted that overestimates of fossil occurrences should lead to underestimates of evolutionary rates. Essentially, as the number of fossil occurrences is increased, the inferred preservation rate will be higher; this increased preservation rate, in turn, will decrease the amount of unobserved evolutionary history that is assumed to have occurred. By reducing the inferred unobserved history through artificial inflation of the preservation rate, we underestimate the number of events that have occurred (speciation or extinction), and the corresponding rate estimates are biased downwards. When the observed number of fossil occurrences is lower than the true number, however, we will tend to overestimate both the unobserved evolutionary history and the corresponding diversification rates.

We also tested the impact of diversification histories that violate the assumptions of the BAMM framework using two scenarios. First, we asked how the recent extinction of a large clade may bias the estimates from BAMM. To simulate this scenario, we chose a single, large, extant clade on each simulated phylogeny and pruned each extant branch back to the time of the most recent fossil occurrence (or dropped it entirely, if it was never observed in the fossil record). For the 2nd scenario, we simulated trees under a diversity-dependent process, and inferred rates from that tree both with and without fossil data. Further details on these scenarios are presented in the Supplementary Materials (Section S8 available on Dryad).

Empirical Analysis

Para efeitos de comparação, incluímos uma análise empírica usando BAMM de uma filogenia de cavalos de 138 pontas recentemente publicada (Cantalapiedra et al. 2017). Escalamos o tempo da árvore de cavalo usando a mesma abordagem de Cantalapiedra et al. (2017), que tende a homogeneizar galhos, assumindo que todos os comprimentos do ramo sejam retirados de uma única distribuição definida por uma taxa de especiação, extinção e preservação (ver Bapst 2013; e os Materiais Suplementares disponíveis na Dryad). Comparamos os resultados do BAMM com os obtidos utilizando-se o PyRate (Silvestro et al. 2014a,b), um quadro de inferência bayesiano que permite a estimativa de parâmetros de diversificação e mudanças de taxa apenas dos dados de ocorrência. Baixamos 1495 ocorrências fósseis (ver Discussão) para as 138 espécies na filogenia do Banco de Dados de Paleobiologia (Alroy et al. 2001) usando a interface R "paleobioDB" (Varela et al. 2015). Simulamos a distribuição posterior das configurações de mudança de taxa usando BAMM com 20.000.000 gerações de amostragem mcmc. A análise pyrate utilizou 1.000.000 gerações de amostragem mcmc.

Resultados

As árvores simuladas variaram de 50 a 4578 pontas (média 770), e a degradação das árvores reduziu o número observado de linhagens para um tamanho médio de 589, 616 e 698 dicas para as menores e mais altas taxas de preservação. Como foram impostas múltiplas histórias de fossilização em cada história evolutiva fixa, o número de ocorrências fósseis é diretamente com a taxa de preservação (números médios de fósseis de 25,4, 255,6 e 2549,9 para as menores e mais altas taxas de preservação). Após 50 milhões de gerações, as cadeias mcmc para pelo menos 95% dos conjuntos de dados em cada esquema de fossilização atingiram nosso critério de convergência alvo, com tamanhos amostrais eficazes para o número de turnos em mais de 200. As simulações do MCMC que não atenderam a esse critério de convergência são, no entanto, incluídas em todas as análises seguintes para evitar viés da subsampleidade, incluindo apenas árvores que apresentaram bom desempenho do MCMC. Os resultados baseados em todas as árvores, mesmo aquelas que não conseguiram alcançar os tamanhos amostrais eficazes desejados, são quase idênticos aos resultados apresentados abaixo e são mostrados no Material Suplementar disponível na Dryad.

=

$=$

Para cada regime de taxa em cada árvore simulada, comparamos as taxas estimadas com as taxas verdadeiras no processo gerador. As taxas estimadas de especiação foram razoavelmente bem correlacionadas com as taxas reais em todos os níveis de preservação. Em todos os regimes de alíquotas e cenários de preservação, a correlação entre as taxas de especiação verdadeiras e estimadas foi de 0,49, embora 70,9% desses regimes contiveram menos de cinco pontas. É improvável que o BAMM tivesse poder suficiente para inferir tais regimes tão pequenos (Rabosky et al. 2017), e as taxas estimadas para a maioria tenderiam a igualar a taxa do processo paterro. Para regimes com cinco ou mais pontas, a correlação aumenta para 0,77 e 0,86 para os somente dependentes e de alta preservação (⁠

ψ=

$\psi =$ 1.0) cenários, com todos os outros cenários de preservação caindo entre esses valores (Fig. 2, coluna esquerda). Para regimes com 20 ou mais pontas, essas correlações sobem para 0,90 e 0,97, respectivamente (Fig. 3, coluna esquerda). As encostas de regressão estão próximas do valor esperado de uma, mas não mostram um aumento consistente com a preservação (somente de subalente: inclinação5+

=

${}_{5+}=$ 0.91 ±

$\pm$ 0.01 s.e. e inclinação20+

${}_{20+}$ 1.01 ±

$\pm$ 0,01 s.e. para regimes com 5 dicas e regimes com 20 dicas, respectivamente; +

$+$+

$+$ψ=

$\psi =$ 0.01: inclinação5+

${}_{5+}$ 0.90 ±

$\pm$ 0.01 s.e. e inclinação20+

${}_{20+}$ 0.99 ±

$\pm$ 0.01 s.e.; ψ=

$\psi =$ 0.1: inclinação5+

${}_{5+}$ 0.86 ±

$\pm$ 0.01 s.e. e inclinação20+

${}_{20+}$ 0.94 ±

$\pm$ 0.01 s.e.; ψ=

$\psi =$ 1: inclinação5+

=

${}_{5+}=$ 0.88 ±

$\pm$ 0.01 s.e. e inclinação20+

=

${}_{20+}=$ 0.96 ±

$\pm$ 0,01 s.e.). As correlações da taxa de especiação melhoram à medida que são considerados regimes com um número mínimo crescente de gorjetas, independentemente da taxa de preservação (Fig. 4).

Figura 2.

Open in new tabDownload slide

BAMM estimates of speciation (⁠

λ)

$\lambda )$⁠, extinction (⁠μ)

$\mu )$⁠, and relative extinction rate (⁠μ

$\mu$/λ)

$\lambda )$ estimates as a function of the true values for rate regimes with at least five tips from trees fossilized with different preservation rates. Each gray point is the estimated rate for a single rate regime. Black points represent mean estimated and inferred rates for all regimes falling within each 2% quantile of true rates. For example, the mean estimated rate of all individual regimes with a true rate between the lowest and the 2% quantile comprise the first black point. Rows give results for extant-only trees (1st row), low (⁠ψ=

$\psi =$ 0.01; 2nd row), moderate (⁠ψ=

$\psi =$ 0.1; 3rd row), and high (⁠ψ=

$\psi =$ 1.0; 4th row). Trees with higher preservation rates, and thus more fossils, have better rate estimates, with the extinction rate showing particular improvement.

$$$$$$$$$$$$$$

Figure 3.

Open in new tabDownload slide

BAMM estimates of speciation (⁠

λ)

$\lambda )$⁠, extinction (⁠μ)

$\mu )$⁠, and relative extinction (⁠μ

$\mu$/λ)

$\lambda )$ rates for rate regimes with at least 20 observed tips from trees fossilized with different preservation rates (⁠ψ)

$\psi )$⁠. Each gray point represents the mean posterior estimate of rates for a regime from the extant-only trees (1st row), low (⁠

ψ

$\psi$=

$=$ 0.01; 2nd row), moderate (⁠ψ=

$\psi =$ 0.1; 3rd row), and high (⁠ψ=

$\psi =$ 1.0; 4th row). Each black point represents the value of all regimes falling within each 2% quantile (e.g., all regimes between the 0th and 2nd percent quantile for the first point). Rates in general are better estimated in the large regimes (as compared to Fig. 2). Increased preservation rates lead to better macroevolutionary rate estimation overall, particularly for extinction.

$$$$$$$$$$$$$$$$$$

Figure 4.

Open in new tabDownload slide

Speciation (⁠

λ)

$\lambda )$ and extinction (⁠μ)

$\mu )$ rate correlations across all regimes as a function of the minimum number of tips present in a particular subset of the data. Each point represents the correlation for speciation (a) or extinction (b) for all rate regimes with the specified minimum number of tips. A minimum regime size (⁠

x

$x$-axis) value of zero indicates that the correlation is computed across all rate regimes, regardless of the number of tips they contain; a minimum value of 10 would restrict this pool of regimes to only those with 10 or more tips. Minimum values of 5 and 20 correspond to results shown in Figs. 2 and 3, respectively. Larger regimes yield better rate estimates for both speciation and extinction. Speciation rates are reliably estimated even in the absence of fossil data, but extinction rate accuracy is markedly improved by the inclusion of fossils.

$$$$$$

Estimates of the extinction rate (⁠

μ)

$\mu )$ are substantially influenced by the preservation rate. Across all regimes and preservation rates, the mean correlation between true and estimated regime-specific rates was 0.14. As for speciation, however, these correlations improve substantially when we exclude the 70.9% of rate regimes that included fewer than five tips. For regimes with five or more tips, correlations ranged from 0.30 for extant-only trees to 0.72 for the ψ=

$\psi =$ 1.0 scenario (Fig. 2). With 20 or more tips, these correlations increase to 0.49 and 0.91, respectively (Fig. 3). As with speciation, the slopes do not vary consistently with preservation, although the error decreases with increasing preservation (extant-only: slopeall

${}_{\mathrm{a}\mathrm{l}\mathrm{l}}$=

$=$ 0.54 ±

$\pm$ 0.02 and slope20+

=

${}_{20+}=$ 1.20 ±

$\pm$ 0.05 for all and regimes with 20 tips; +

$+$ψ=

$\psi =$ 0.01: slopeall

=

${}_{\mathrm{a}\mathrm{l}\mathrm{l}}=$ 0.47 ±

$\pm$ 0.02 and slope20+

=

${}_{20+}=$ 0.95 ±

$\pm$ 0.05; ψ=

$\psi =$ 0.1: slopeall

=

${}_{\mathrm{a}\mathrm{l}\mathrm{l}}=$ 0.38 ±

$\pm$ 0.01 and slope20+

=

${}_{20+}=$ 0.78 ±

$\pm$ 0.02; ψ=

$\psi =$ 1: slopeall

=

${}_{\mathrm{a}\mathrm{l}\mathrm{l}}=$ 0.50 ±

$\pm$ 0.01 and slope20+

=

${}_{20+}=$ 0.95 ±

$\pm$ 0.01).

The relative extinction rates (⁠

μ

$\mu$/λ)

$\lambda )$ also improved with increasing preservation, with the estimates from extant-only simulations only weakly correlated with the true relative extinction level (⁠r=

$r=$ 0.41 for 5 tips, 0.53 for 20 tips) and again increasing with increased preservation (+

$+$+

$+$⁠r=

$r=$ 0.67 for regimes with 5 tips; 0.87 for regimes with 20 tips; Figs. 2 and 3; right column). As for speciation, extinction rate correlations improve as we consider subsets of regimes with increasing numbers of tips, but accuracy of rate estimates is heavily influenced by the preservation rate (Fig. 4b). Preservation rates were estimated with high accuracy, with tight distributions around the generating values at +

$+$>

$>$ψ=

$\psi =$ 0.01, 0.1, and 1.0 (the 5 and 95% quantiles of the posterior estimate for ψ=

$\psi =$ 0.01 are 

ψ

5

=

${\psi }_5=$0.007 and 

ψ

95

=

${\psi }_{95}=$ 0.015; 

ψ

5

=

${\psi }_5=$ 0.09 and 

ψ

95

=

${\psi }_{95}=$ 0.12 for ψ=

$\psi =$ 0.1; 

ψ

5

=

${\psi }_5=$ 0.96 and

ψ

95

=

${\psi }_{95}=$ 1.04 for ψ=

$\psi =$ 1; Fig. 5).

Figure 5.

Open in new tabDownload slide

Marginal posterior distributions for the preservation rate (⁠

ψ)

$\psi )$ for trees simulated with preservation rates of 0.01, 0.1, and 1.0. Black arrows denote the true value of 

ψ

$\psi$⁠, while gray bars show the distribution of mean 

ψ

$\psi$ across all 200 trees for each value.

$$$$$$

As has been shown previously (Rabosky 2014), BAMM tends to be relatively conservative in the number of regimes on the tree even with the high prior on the number of shifts used in our analyses (Fig. 6). Using stepwise model selection with a Bayes factor threshold of 20, we identified spurious shifts in fewer than 0.5% of trees (Fig. 6). We tabulated the marginal posterior probability of a rate shift for each branch on which a (true) rate shift occurred. These marginal probabilities varied substantially but generally increased with the number of descendent tips across all preservation regimes (Fig. 7a). The marginal probability of a shift along a particular branch increases with both the magnitude of the rate increase and with the number of descendant lineages (Fig. 7b). Thus, BAMM substantially underestimates the number of rate shifts (Fig. 6), but shifts that are not detected tend to be those leading to small clades and/or those that represent small changes in evolutionary rates (Fig. 7). Tip partitioning accuracy in BAMM v2.6 was high overall, and did not vary systematically with preservation rate (Fig. 7c).

Figure 6.

Open in new tabDownload slide

The true number of shifts in the degraded trees plotted against the mean number of inferred shifts (number of non-root rate regimes best supported by stepwise model selection using Bayes factors) for preservation rates (⁠

ψ)

$\psi )$ of 0.01, 0.1, and 1.0. Each point represents one of the 200 trees. The identity line is shown for reference, and the estimates slopes are very low (0.08, 0.085, and 0.18). The correlation between true and estimated numbers of shifts are 0.65, 0.73, and 0.8 for the preservations rates from lowest to highest. Plots in the top row show the results on the identical axes, while the plots in the 2nd row show the same data with the y

$y$-axis rescaled to better show the correlation between the true and estimated number of shifts. Most simulated shifts are small in terms of the number of lineages and they may also reflect a relatively minor change in rates (e.g., a shift from speciation of 0.1–0.11). As such, a method focused on summarizing the major regimes of the tree (i.e., those with strong support) is expected to underestimate the number of shifts. This is especially true when a stringent Bayes factor threshold is used (see SOM for the mean number of shifts from the posterior distribution).

$$$$

Figure 7.

Open in new tabDownload slide

Topological accuracy of inferred shift locations. a) Marginal posterior probability of each rate shift as a function of the number of tips in the shift regime. A value of 0.50 would imply that 50% of samples from the posterior recovered a shift in the same topological location (e.g., same branch) as the true shift. In general, we expect the marginal shift probability to be somewhat less than 1.0 even for well-supported shifts, because the evidence for a given shift can be smeared across several adjacent branches (Rabosky 2014). Red horizontal lines denote median marginal probabilities of rate shifts within sliding windows of arbitrary size (see text for details). b) Marginal probability of each rate shift as a function of both the proportional change in speciation rate (⁠

λNEW

${\lambda }_{\mathrm{N}\mathrm{E}\mathrm{W}}$/λOLD

)

${\lambda }_{\mathrm{O}\mathrm{L}\mathrm{D}})$ as well as the number of descendant tips in the regime. Darker points represent higher marginal probabilities. Shift locations are inferred with the highest accuracy is highest when large rate increases occur and/or when shifts lead to large numbers of descendant tips (c) Accuracy with which BAMM infers cohort membership, computed as the mean probability of assigning a given taxon pair to the correct grouping relationship (e.g., “same shift regime”, “different shift regime”). Red points represent the mean cohort accuracy of BAMM analyses across four preservation rates as a function of the number of shifts. Black squares represent the corresponding cohort assignment accuracy expected under random shift placement; see text for details.

${}_{\mathrm{}\mathrm{}\mathrm{}}$${}_{\mathrm{}\mathrm{}\mathrm{}}$

The apparent conservatism of BAMM with respect to the true number of shifts need not be a problem: we have previously shown that, under the Poisson process, the vast majority of rate shifts lead to small rate regimes (⁠

<

$<$5 tips), and that shift regimes of this size are unlikely to be inferred by any method (Rabosky et al. 2017). Because the regimes are small (e.g., contain minimal data), their likelihood is similar across a broad range of parameter values, and there is insufficient information in the data to infer their existence. As pointed out by one reviewer, our results also raise a paradox: how can the BAMM-estimated rates perform well, despite the massive underestimate of the number of shifts? The solution to this paradox lies in the fact that small rate shifts are missed by BAMM, but dominant rate classes across the tree are captured by BAMM with reasonable accuracy. Hence, BAMM might correctly infer only 3 of 15 rate regimes in a given phylogeny, but those three inferred regimes can govern the majority of branches in the tree (e.g., 12 shifts leading to one or two taxa, but the other three shifts all leading to 50 or more taxa).

Sensitivity of the Method to Model Violations

We found that the implementation performed well even when assumptions of homogeneous fossilization rates through time were violated. After imposing temporally-varying preservation rates on the data, we found that diversification rates were still estimated with high accuracy (Fig. 8). The estimates of both speciation and relative extinction (Fig. 8a,b) were reasonably correlated with the generating values. We also found very few instances where the number of inferred shifts was higher than the generating model (Fig. 8c). Consistent with the time-homogenous preservation simulations, 78% of the tips were accurately partitioned among the rate regimes (Fig. 8d).

Figure 8.

Open in new tabDownload slide

Performance of BAMM when the assumption of time-homogeneous fossil preservation is violated. Speciation (⁠

λ

$\lambda$⁠; a), relative extinction (⁠μ

$\mu$/λ

$\lambda$⁠; b), inferred number of shifts (c) and cohort accuracy (d) are shown for trees fossilized using a stage-specific simulation and analyzed using the constant-preservation rate assumption of BAMM; results illustrate regimes with five or more tips, for comparison to Fig. 2.

$$$$$$

Underestimating the number of fossil occurrences inflates the speciation and extinction rates, while overestimating the number of occurrences deflates them (Fig. 9). Individual estimates of 

λ

$\lambda$ may be off by as much as 50% when misspecification is exceptionally high (Fig. 9a), while μ

$\mu$ may vary by substantially more (Fig. 9b). However, in aggregate across all 100 trees, the median error at each level of occurrence misspecification is fairly low. Encouragingly, the correlations of tree-wide rate estimates for these constant rate trees are high even when the fossil occurrence numbers are off by a factor of ten (correlation for speciation: 0.95 to 0.89; correlation for extinction: 0.93 to 0.89). This suggests that, if the number of fossil occurrences is known only with large uncertainty, the relative rates across the phylogeny are likely to be more accurate than the absolute values of the rates (e.g., shifts and temporal heterogeneity may still be detected reliably, even if the rate magnitude is biased).

Figure 9.

Open in new tabDownload slide

Bias in the estimation of macroevolutionary rates when the number of fossil occurrences was misspecified. The minimum number of fossil occurrences (F

MIN

)

${}_{\mathrm{M}\mathrm{I}\mathrm{N}})$ is equal to the number of observed extinct tips, and the gap between the minimum and true number of fossil occurrences (FGAP

)

${}_{\mathrm{G}\mathrm{A}\mathrm{P}})$ was used to analyzed different trees at comparable levels of misspecification. The bias in speciation (a) and extinction (b) show a similar pattern of inflated rate estimates when fewer than the true number of observed fossils are specified, and deflated estimates when the number of fossils were exaggerated. Bold points represent the median bias, and red lines indicated the true number of fossils (vertical) and no bias in rate estimation (horizontal). The degrees of fossil misspecification are: 1) FMIN

${}_{\mathrm{M}\mathrm{I}\mathrm{N}}$⁠, 2) FMIN

+

${}_{\mathrm{M}\mathrm{I}\mathrm{N}}+$ 0.25 ×

$\times$ FGAP

${}_{\mathrm{G}\mathrm{A}\mathrm{P}}$⁠, 3) FMIN

+

${}_{\mathrm{M}\mathrm{I}\mathrm{N}}+$ 0.5 ×

$\times$ FGAP

${}_{\mathrm{G}\mathrm{A}\mathrm{P}}$⁠, 4) FTRUE

${}_{\mathrm{T}\mathrm{R}\mathrm{U}\mathrm{E}}$⁠, 5) FTRUE

+

${}_{\mathrm{T}\mathrm{R}\mathrm{U}\mathrm{E}}+$ FGAP

${}_{\mathrm{G}\mathrm{A}\mathrm{P}}$⁠, and 6) 10 ×

$\times$ FTRUE

${}_{\mathrm{T}\mathrm{R}\mathrm{U}\mathrm{E}}$⁠. Trees were simulated using constant rates through time, and the inferred values are for the entire tree.

${}_{\mathrm{}\mathrm{}\mathrm{}}$${}_{\mathrm{}\mathrm{}\mathrm{}}$${}_{\mathrm{}\mathrm{}\mathrm{}}$${}_{\mathrm{}\mathrm{}\mathrm{}}$$$${}_{\mathrm{}\mathrm{}\mathrm{}}$${}_{\mathrm{}\mathrm{}\mathrm{}}$$$${}_{\mathrm{}\mathrm{}\mathrm{}}$${}_{\mathrm{}\mathrm{}\mathrm{}\mathrm{}}$${}_{\mathrm{}\mathrm{}\mathrm{}\mathrm{}}$${}_{\mathrm{}\mathrm{}\mathrm{}}$$$${}_{\mathrm{}\mathrm{}\mathrm{}\mathrm{}}$

We explored BAMM’s performance under two scenarios where the diversification process in the simulation model violated the assumptions of the inference model. For the diversity-dependent simulations, BAMM performed much better with the inclusion of fossil information (Supplementary Fig. S11 available on Dryad): both speciation and extinction rates were estimated with much greater accuracy when fossils were included, relative to the extant-only analysis. For the mass extinction scenario, we observed consistent biases in rate estimates when fossil and extant data were combined, and relative extinction rates were consistently biased upwards (Supplementary Fig. S10 available on Dryad).

Unobserved shifts seem to have a negligible impact on rate estimation. Neither the slope nor correlation coefficient of the extinction and speciation rates for branches within trees varies with the proportion of unobserved rate shifts (Fig. 10), nor with the absolute number of unobserved shifts (see Supplementary Material available on Dryad). In contrast to the regime-specific assessments presented above, within-tree measures are difficult to interpret under rate-shift models like BAMM. If a method does not infer rate heterogeneity within the tree, the metrics (especially the slope) are expected to perform poorly. Rabosky et al. (2017) discusses this issue in detail, documenting how the metrics may be misleading even when rates on almost every branch are estimated with high accuracy. Regardless, there is no obvious trend between the proportion of unobserved shifts within a tree and the overall metrics of performance (Fig. 10).

Figure 10.

Open in new tabDownload slide

Metrics of within-tree rate accuracy as a function of the proportion of shifts in the tree that are unobserved, for trees with at least 50 observed tips. For the branches in each tree we evaluated the inferred and generating rate of speciation and extinction, calculated the slope (1st and 3rd rows), and Pearson correlation (2nd and 4th rows) between them, and show them relative to the proportion of rate shifts in the tree that have no sampled taxa (see Fig. 1). The relative density of values across all preservation levels is shown to the right. Both metrics for both rates have their highest density across all trees and rate regimes near the optimal value of 1 (mean slopes of 0.67 and 0.41 for speciation and extinction; mean correlations of 0.77 and 0.48 for speciation and extinction). Neither rate shows a relationship between the within-tree metric and the proportion of unobserved shifts suggesting that unobserved shifts do not bias rate estimates in these simulations. Filled circles denote trees where BAMM identified substantial evidence for rate variation (Bayes factor evidence 20 for one or more shifts); open circles are trees that lacked evidence for rate variation. Note that branch-specific slopes and correlations as presented above are expected to equal zero if BAMM fails to detect diversification rate variation (see Rabosky et al. 2017, Fig. 10). For additional visualizations demonstrating the lack of a relationship between performance and unobserved shifts, see the Supplementary Materials available on Dryad.

>

$>$

$$

Empirical Analysis

Our empirical analysis of the horse phylogeny yielded evidence for a major shift in rates leading up to the extant genus Equus, and also that speciation and extinction rates are approximately equal across the tree and through time (Fig. 11). Net diversification was found to be near-zero across the horse phylogeny (mean net div. 

=−

$=-$0.002; see Fig. 11), but lineage turnover rates were substantial (mean turnover 0.7). This result is concordant with long-standing hypotheses about the history of life, where speciation is only slightly elevated relative to extinction, and only in some clades (e.g., Raup 1994). However, it is worth noting that the BAMM estimates assume that rates have been constant through time within rate regimes, and it is unlikely that this assumption is strictly true in practice. It is possible, for example, that the near-zero net diversification rates estimated for equids represent the long-term average of a clade that experienced a pronounced phase of positive net diversification followed by negative net diversification rates that led to extinction of many lineages (Quental and Marshall 2013; Cantalapiedra et al. 2017). BAMM’s assumption of time-homogeneous diversification means, in practice, that temporal rate variation will be smeared across the duration of a particular rate regime. We caution users of BAMM from overinterpreting rates at specific points in time if time-constancy of rates is assumed within rate regimes.=

$=$

Figure 11.

Open in new tabDownload slide

Mean net diversification rates per branch (left), and mean speciation and extinction rates through time from BAMM and PyRate (right) for the horse phylogeny. PyRate and BAMM estimates disagree with respect to the overall turnover and net diversification rates through time, although estimates of speciation are in accordance at the beginning of the horse phylogeny (⁠

∼

$\sim$20–15 Ma) and extinction estimates are comparable toward the present (⁠∼

$\sim$3–0 Ma).

$$$$

PyRate likewise infers an increase in extinction through time, but differs from the BAMM results in several ways. Extinction as estimated with PyRate was substantially lower than speciation during the earliest portion of horse evolution, but speciation dropped substantially at roughly 15 Ma and rebounded slightly at 6 Ma. The estimates of speciation and extinction from PyRate differ substantially from those estimated by BAMM, with greater volatility in the PyRate estimates of net diversification; BAMM inferred a more consistent near-zero level of net diversification. Notably, the turnover rates estimated by BAMM are much higher than those from PyRate. Both methods estimate high preservation rates, although as speciation, extinction, and preservation are estimated jointly, the two estimates are different, with BAMM estimating a higher rate (⁠

ψ=

$\psi =$ 3.69 ±

$\pm$ 0.10) than PyRate (⁠ψ=

$\psi =$ 2.27 ±

$\pm$ 0.07).

Discussion

We have described an extension of the BAMM framework for modeling macroevolutionary dynamics that leverages recent advances in modeling the fossilized birth–death process (Stadler 2010; Ronquist et al. 2012; Heath et al. 2014; Didier et al. 2012, 2017) to estimate rates of speciation, extinction, and preservation from phylogenies that include observations of fossil (non-extant) lineages. We found that BAMM was able to infer diversification rates, preservation rates, and the number and location of rate shifts with reasonable accuracy across a range of fossilization scenarios. The inference model implemented in BAMM is robust to at least some violations of its assumptions, although we necessarily explored a limited subset of many possible ways that model assumptions can be violated. Given debates over the inference of extinction from molecular phylogenies (Rabosky 2010, 2016; Beaulieu and O’Meara 2015), it is noteworthy that inclusion of even a small number of fossils (⁠

<

$<$10% of the total tips) substantially improves the estimation of both the absolute and relative extinction rates in comparison to analyses of the same data restricted to extant tips only. BAMM’s likelihood is conditioned on the absence of unobserved rate shifts, an assumption that has been criticized in the recent literature (Moore et al. 2016). However, we find no evidence that these unknown quantities affect inference (Fig. 10). Our results are thus consistent with Rabosky et al.’s (2017) demonstration that unobserved rate shifts are unimportant for empirically relevant parameterizations of the diversification process.

Given the preservational biases inherent to the fossil record, our results showing that BAMM is at least somewhat robust to temporal heterogeneity in preservation rates and uncertainty in the number of fossil occurrences are encouraging (Figs. 8 and 9). Changes in the preservation rate through time are likely to be more extreme in real data sets than considered by our study. However, as long as the average preservation through time is comparable to values used in our simulations, we expect that the inclusion of fossil data will improve estimates of diversification parameters. Incorrect specification of the number of fossil occurrences biases the magnitude of rate estimates, but the true and estimated rates remain highly correlated even with substantial misspecification of the true number of fossils. This suggests that relative rates within a phylogeny may largely be robust to misspecification of occurrence counts, even if the absolute values of those rates are biased. Nonetheless, assumptions about the nature of fossil preservation and lineage comparability between neontological and paleontological datatypes should be treated with caution in future empirical work.

In the case of clade-wide mass extinction (Supplementary Fig. S10 available on Dryad), the bias in BAMM’s extinction estimates is easily understood. BAMM has no ability to accommodate mass extinctions directly, but it is nonetheless forced to model a situation where a large number of taxa became extinct at a particular point in time. Hence, the background extinction rate for the focal clade must increase to accommodate the extinction of many taxa, even though the clade may have had low extinction throughout much of its history. In effect, BAMM smears the effects of the mass extinction across earlier periods of time when background extinction was lower.

Because the current implementation of the model does not accommodate time-varying diversification rates, it is possible that resulting time-smeared rate estimates will be an inadequate description of the true diversification process. Similarly, for any clades that have gone extinct in their entirety, net diversification rates will generally be estimated as near zero throughout the history of the clade. This potential for time-averaging of rates is a substantial concern for interpretation that must be considered carefully by users of the method, at least until extensions have been developed that can better accommodate time-varying diversification processes.

Sampling Biases

Paleontologists have long recognized that apparent long-term trends in macroevolutionary dynamics can result from biases in the rock record, and in particular the temporal variation in the total amount of sedimentary rock (Raup 1972, 1976; Sepkoski et al. 1981; Foote and Raup 1996; Peters and Foote 2001; Peters 2006; Holland 2016). Even if total species richness is constant through time, intervals with greater amounts of exposed sedimentary rock volume, or with the same amount of rock distributed across a greater geographic extent, will generally appear to be more diverse (Raup 1972). A general trend of increasing preservation potential through time has long been documented (e.g., Raup 1976), and if this trend were universally monotonic it would be relatively easy to model. However, preservation rates have varied geographically throughout Earth’s history (Peters 2005). At any given time in the past, most of the habitable environments that existed are unpreserved. This facet of the rock record becomes particularly problematic when a clade endemic to a region either appears or disappears due to a shift in preservation potential. If 

ψ

$\psi$ does not vary to accommodate the change in local sampling rates, then the sudden change in diversity may be incorrectly inferred to reflect a change in the diversification parameters.

Within-lineage preservation heterogeneity can also impact macroevolutionary inference, and has been demonstrated to influence estimates of diversification (Liow et al. 2010b; Silvestro et al. 2014a,b). Factors such as geographic occupancy (Foote et al. 2007, 2016) and expected abundance (Stanley 1979) are thought to vary with the age of a lineage, and both can influence preservation potential. The probability of preservation may also vary among lineages within a phylogeny. Many factors, from body size and shell composition to habitat preferences and foraging strategies, influence the preservation potential of a taxon (Valentine 1989; Behrensmeyer and Hook 1992; Behrensmeyer et al. 2000; Valentine et al. 2006). Furthermore, the magnitude and directionality of these biases can vary with time and underlying sedimentology (Kidwell 2005; Foote et al. 2015). The BAMM extension described in this article incorporates a very simple model of preservation, such that users merely need to specify the total number of fossil occurrences associated with the lineages that are included in the phylogeny. Accommodating changes in preservation potential through time, among lineages, and across geographic regions is an obvious future extension for models like BAMM that integrate extant and fossil data.

Comparability of Neontological and Paleontological Data

Many challenges remain in linking paleontological and neontological data, as the very concept of a lineage differs between the two datatypes. Fossil taxa are identified based on preserved morphological traits, such that that a single fossil morphospecies may be equivalent to multiple biological species that are distinguishable by multiple traits that are unlikely to fossilize (e.g., behavior, coloration, soft-tissue anatomy, vocalizations). Further, it is difficult in the fossil record to discriminate between a single lineage sampled repeatedly through time (a series of anagenetic morphospecies) and multiple, closely-related but separate lineages present in different stratigraphic horizons. Hence, inferences that use method described in this article are subject to an additional source of uncertainty that reflects different species concepts (either in theory or in practice) between neontological and paleontological data.

How Many Fossil Occurrences?

The model we have implemented assumes that fossil occurrences are generated under a simple birth–death-fossilization process, which implies that occurrences represent the sampling of a specific lineage at a specific point in time. We suggest that the number of occurrences should reflect the number of temporally-unique occurrences of identifiable lineages, or the number of stratigraphically unique occurrences of both extinct and extant species-level taxa that are included in the phylogeny. The number of temporally-unique identifiable occurrences differs from raw counts of the total number of fossil occurrences of the clade in several important ways. In BAMM, taxonomically indeterminate records belonging to a clade (e.g., “Aves indet.”) should not be used, as there is no way to know whether the fossil corresponds to an observed lineage within the focal clade (either an extinct or extant form represented by a lineage on the phylogeny), or if it belongs to an as-yet unrecognized branch within the clade (e.g., a new species not yet included in the phylogeny). Assuming we are generally working with species-level phylogenies, one should also avoid using occurrences that have not been resolved below the genus level in BAMM, as they might represent lineages that are not present in the tree. Furthermore, multiple individuals of the same species collected from a single bone or shell bed should be treated as a single occurrence.

As an example of how researchers might specify an appropriate number of fossil occurrences for a BAMM analysis, we will consider how occurrences might be tabulated for a recent phylogenetic study of extant and extinct canid mammals (Matzke and Wright 2016). A simple search (17 February 2017) on “Canidae” from the Paleobiology Database (Alroy et al. 2001; Varela et al. 2015) gave at least 10,660 fossil observations (unique records). Of this set, 2946 represented stratigraphically unique occurrences, and 2107 of those observations were resolved to the species level. The data set analyzed by Matzke and Wright (2016) contains 142 species-level taxa, and we identified 1837 stratigraphically-unique occurrences of those taxa in the database. For the simple (time-homogeneous) preservation model described in this article, we would specify 1837 occurrences for a BAMM analysis of the canid phylogeny. The Paleobiology Database allows rapid searching for occurrences that match only the set of taxa present in a given phylogeny, and querying the output for duplicate localities can be done quickly (see Dryad scripts for our horse example).

At present, the model assumes that fossil occurrences are comparable through time and across clades. This assumption is clearly violated in practice, particularly when we consider large clades whose constituent species differ in key biological traits that influence their probability of preservation. To illustrate this point, consider isolated vertebrate teeth, a common type of fossil. Crocodilians shed their teeth continuously through life, although their teeth are rarely identifiable to the species level. On the other hand, mammals only have two sets of teeth in their life, so they produce fewer fossils, but each tooth is far more likely to be mapped to an observed lineage on a phylogeny. These differences make it difficult to equate the occurrence of a single crocodilian tooth to the occurrence of a mammalian molar. An obvious extension to the modeling framework presented here will be to incorporate lineage-specific variation in preservation rates, which will help account for such clade-specific differences in preservation potential.

Conclusions

The fossilized birth–death process has already shown great promise for inferring and time-scaling phylogenies using all available evidence. Recent advances in extending the model to incorporate piece-wise time variation (Gavryushkina et al. 2014, 2017) and non-uniform incomplete sampling (Zhang et al. 2016) have greatly improved the ability to estimate phylogenies that include both extinct and extant lineages. Herein, we have shown that BAMM performs well on dated phylogenies that include a mix of extant and extinct tips (or extinct-only tips), providing a useful tool for characterizing diversification rate variation across the expanding pool of time-calibrated phylogenies that include extinct lineages. This approach will help researchers test hypotheses about macroevolutionary rate control in clades that may lack sufficient fossil data for non-phylogenetic (occurrence-only) paleobiological inferences.

The results of this study and others (e.g., Didier et al. 2012) demonstrate that the inclusion of fossil data improves the estimation of both speciation and extinction rates from phylogenies. The assumptions of the model we describe are simplistic, and numerous challenges remain to be addressed. However, we have found that results are reasonably robust to several key violations of model assumptions. The method described here is a first step toward integrating information from paleontological and neontological data sets to better understand the extent to which the dynamics of species diversification have varied through time and among lineages.

Supplementary Material

Data available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.50m70.

Funding

This work was supported by the U.S. National Science Foundation [NSF-DEB-1256330 to D.L.R.], by a fellowship from the David and Lucile Packard Foundation [D.L.R.], and by a VICI grant [R.S.E.] from the Netherlands Organization for Scientific Research (NWO).

Acknowledgments

We thank R. Ree and three anonymous reviewers for comments on the article.

References

Alfaro

 

M.E.

, 

Santini

 

F.

, 

Brock

 

C.

, 

Alamillo

 

H.

, 

Dornburg

 

A.

, 

Rabosky

 

D.L.

, 

Carnevale

 

G.

, 

Harmon

 

L.J.

 

2009

. 

Nine exceptional radiations plus high turnover explain species diversity in jawed vertebrates.

 

Proc. Natl. Acad. Sci. USA.

 

106

:

13410

–

13414

.

Google Scholar

Crossref

 

Alroy

 

J.

 

1999

. 

The fossil record of North American mammals: evidence for a paleocene evolutionary radiation.

 

Syst. Biol.

 

48

:

107

–

118

.

Google Scholar

Crossref

PubMed

 

Alroy

 

J.

 

2008

. 

Dynamics of origination and extinction in the marine fossil record.

 

Proc. Natl. Acad. Sci. USA.

 

105

:

11536

–

11542

.

Google Scholar

Crossref

 

Alroy,

 

J.

 

2010

. 

The shifting balance of diversity among major marine animal groups.

 

Science.

 

329

: 

1191

–

1194

.

Google Scholar

Crossref

PubMed

 

Alroy

 

J.

, 

Marshall

 

C.R.

, 

Bambach

 

R.K.

, 

Bezusko

 

K.

, 

Foote

 

M.

, 

Fürsich

 

F.T.

, 

Hansen

 

T.A.

, 

Holland

 

S.M.

, 

Ivany

 

L.C.

, 

Jablonski

 

D.

, 

Jacobs

 

D.K.

, 

Jones

 

D.C.

, 

Kosnik

 

M.A.

, 

Lidgard

 

S.

, 

Low

 

S.

, 

Miller

 

A.I.

, 

Novack-Gottshall

 

P.M.

, 

Olszewski

 

T.D.

, 

Patzkowsky

 

M.E.

, 

Raup

 

D.M.

, 

Roy

 

K.

, 

Sepkoski

 

J.J.

, 

Sommers

 

M.G.

, 

Wagner

 

P.J.

, 

Webber

 

A.

 

2001

. 

Effects of sampling standardization on estimates of Phanerozoic marine diversification.

 

Proc. Natl. Acad. Sci. USA.

 

98

:

6261

–

6266

.

Google Scholar

Crossref

 

Bapst,

 

DW.

 

2013

. 

A stochastic rate-calibrated method for time-scaling phylogenies of fossil taxa.

 

Methods Ecol. Evol.

 

4

: 

724

–

733

.

Google Scholar

Crossref

 

Beaulieu

 

J.M.

, 

O’Meara

 

B.C.

 

2015

. 

Extinction can be estimated from moderately sized molecular phylogenies.

 

Evolution.

 

69

:

1036

–

1043

.

Google Scholar

Crossref

PubMed

 

Behrensmeyer,

 

A.K.

, 

Hook,

 

R.W.

 

1992

. 

Paleoenvironmental contexts and taphonomic modes.

 In: 

Behrensmeyer,

 

A.K.,

 

Damuth,

 

J.D.,

 

DiMichele,

 

W.A.,

 

Potts,

 

R.,

 

Sues,

 

H.D.,

 

Wing,

 

S.L.,

 editors. 

Terrestrial ecosystems through time

. 

Chicago

: 

University of Chicago Press

. Pp. 

15

–

136

.

Google Scholar

 

Behrensmeyer

 

A.K.

, 

Kidwell

 

S.M.

, 

Gastaldo

 

R.A.

 

2000

. 

Taphonomy and paleobiology.

 

Paleobiology.

 

26

:

103

–

147

.

Google Scholar

Crossref

 

Benton

 

M.J.

, 

Pearson

 

P.N.

 

2001

. 

Speciation in the fossil record.

 

Trends Ecol. Evol.

 

16

:

405

–

411

.

Google Scholar

Crossref

PubMed

 

Cantalapiedra

 

J.L.

, 

Prado

 

J.L.

, 

Fernández

 

M.H.

, 

Alberdi

 

M.T.

 

2017

. 

Decoupled ecomorphological evolution and diversification in Neogene-Quaternary horses.

 

Science.

 

355

:

627

–

630

.

Google Scholar

Crossref

PubMed

 

Didier

 

G.

, 

Fau

 

M.

, 

Laurin

 

M.

 

2017

. 

Likelihood of tree topologies with fossils and diversification rate estimation.

 

Syst. Biol.

 

66

:

964

–

987

.

Google Scholar

Crossref

PubMed

 

Didier

 

G.

, 

Royer-Carenzi

 

M.

, 

Laurin

 

M.

 

2012

. 

The reconstructed evolutionary process with the fossil record.

 

J. Theor. Biol.

 

315

:

26

–

37

.

Google Scholar

Crossref

PubMed

 

Etienne

 

R.S.

, 

Haegeman

 

B.

 

2012

. 

A conceptual and statistical framework for adaptive radiations with a key role for diversity dependence.

 

Am. Nat.

 

180

:

E75

–

E89

.

Google Scholar

Crossref

PubMed

 

Etienne

 

R.S.

, 

Haegeman

 

B.

, 

Stadler

 

T.

, 

Aze

 

T.

, 

Pearson

 

P.N.

, 

Purvis

 

A.

, 

Phillimore

 

A.B.

 

2012

. 

Diversity-dependence brings molecular phylogenies closer to agreement with the fossil record.

 

Proc. R. Soc. Lond. B Biol. Sci.

 

279

:

1300

–

1309

.

Google Scholar

Crossref

 

FitzJohn

 

R.

 

2010

. 

Quantitative traits and diversification.

 

Syst. Biol.

 

59

:

619

–

633

.

Google Scholar

Crossref

PubMed

 

FitzJohn

 

R.G.

, 

Maddison

 

W.P.

, 

Otto

 

S.P.

 

2009

. 

Estimating trait-dependent speciation and extinction rates from incompletely resolved phylogenies.

 

Syst. Biol.

 

58

:

595

–

611

.

Google Scholar

Crossref

PubMed

 

Foote

 

M.

, 

Crampton

 

J.S.

, 

Beu

 

A.G.

, 

Marshall

 

B.A.

, 

Cooper

 

R.A.

, 

Maxwell

 

P.A.

, 

Matcham

 

I.

 

2007

. 

Rise and fall of species occupancy in Cenozoic fossil mollusks.

 

Science.

 

318

:

1131

–

1134

.

Google Scholar

Crossref

PubMed

 

Foote

 

M.

, 

Crampton

 

J.S.

, 

Beu

 

A.G.

, 

Nelson

 

C.S.

 

2015

. 

Aragonite bias, and lack of bias, in the fossil record: lithological, environmental, and ecological controls.

 

Paleobiology.

 

41

:

245

–

265

.

Google Scholar

Crossref

 

Foote

 

M.

, 

Raup

 

D.M.

 

1996

. 

Fossil preservation and the stratigraphic ranges of taxa.

 

Paleobiology.

 

22

:

121

–

140

.

Google Scholar

Crossref

PubMed

 

Foote

 

M.

, 

Ritterbush

 

K.A.

, 

Miller

 

A.I.

 

2016

. 

Geographic ranges of genera and their constituent species: structure, evolutionary dynamics, and extinction resistance.

 

Paleobiology.

 

42

:

269

–

288

.

Google Scholar

Crossref

 

Gavryushkina

 

A.

, 

Heath

 

T.A.

, 

Ksepka

 

D.T.

, 

Stadler

 

T.

, 

Welch

 

D.

, 

Drummond

 

A.J.

 

2017

. 

Bayesian total-evidence dating reveals the recent crown radiation of penguins.

 

Syst. Biol.

 

66

:

57

–

73

.

Google Scholar

PubMed

 

Gavryushkina

 

A.

, 

Welch

 

D.

, 

Stadler

 

T.

, 

Drummond

 

A.J.

 

2014

. 

Bayesian inference of sampled ancestor trees for epidemiology and fossil calibration.

 

PLoS Comput. Biol.

 

10

:

e1003919

.

Google Scholar

Crossref

PubMed

 

Heath

 

T.A.

, 

Huelsenbeck

 

J.P.

, 

Stadler

 

T.

 

2014

. 

The fossilized birth–death process for coherent calibration of divergence-time estimates.

 

Proc. Natl. Acad. Sci. USA.

 

111

:

E2957

–

E2966

.

Google Scholar

Crossref

 

Holland

 

S.M.

 

2016

. 

The non-uniformity of fossil preservation.

 

Philos. Trans. R. Soc. Lond. B Biol. Sci.

 

371

:

20150130

.

Google Scholar

Crossref

PubMed

 

Hunt

 

G.

, 

Slater

 

G.

 

2016

. 

Integrating paleontological and phylogenetic approaches to macroevolution.

 

Annu. Rev. Ecol. Evol. Syst.

 

47

:

189

–

213

.

Google Scholar

Crossref

 

Kidwell

 

S.M.

 

2005

. 

Shell composition has no net impact on large-scale evolutionary patterns in mollusks.

 

Science.

 

307

:

914

–

917

.

Google Scholar

Crossref

PubMed

 

Ksepka

 

D.T.

, 

Parham

 

J.F.

, 

Allman

 

J.F.

, 

Benton

 

M.J.

, 

Carrano

 

M.T.

, 

Cranston

 

K.A.

, 

Donoghue

 

P.C.J.

, 

Head

 

J.J.

, 

Hermsen

 

E.J.

, 

Irmis

 

R.B.

, 

Joyce

 

W.G.

, 

Kohli

 

M.

, 

Lamm

 

K.D.

, 

Leehr

 

D.

, 

Patané

 

J.L.

, 

Polly

 

P.D.

, 

Phillips

 

M.J.

, 

Smith

 

N.A.

, 

Smith

 

N.D.

, 

Tuinen

 

M.V.

, 

Ware

 

J.L.

, 

Warnock

 

R.C.M.

 

2015

. 

The Fossil Calibration Database—a new resource for divergence dating.

 

Syst. Biol.

 

64

:

853

–

859

.

Google Scholar

Crossref

PubMed

 

Liow

 

L.H.

, 

Quental

 

T.B.

, 

Marshall

 

C.R.

 

2010a

. 

When can decreasing diversification rates be detected with molecular phylogenies and the fossil record? Syst.

 

Biol.

 

59

:

646

–

659

.

Google Scholar

 

Liow

 

L.H.

, 

Skaug

 

H.J.

, 

Ergon

 

T.

, 

Schweder

 

T.

 

2010b

. 

Global occurrence trajectories of microfossils: environmental volatility and the rise and fall of individual species.

 

Paleobiology.

 

36

:

224

–

252

.

Google Scholar

Crossref

 

Maddison

 

W.P.

, 

Midford

 

P.E.

, 

Otto

 

S.P.

 

2007

. 

Estimating a binary character’s effect on speciation and extinction.

 

Syst. Biol.

 

56

:

701

–

710

.

Google Scholar

Crossref

PubMed

 

Matzke,

 

N. J.,

 

Wright,

 

A.

 

2016

. 

Inferring node dates from tip dates in fossil Canidae: the importance of tree priors.

 

Biol. Lett.

 

12

:

20160328

Google Scholar

Crossref

PubMed

 

Mitchell

 

J.S.

, 

Makovicky

 

P.J.

 

2014

. 

Low ecological disparity in Early Cretaceous birds.

 

Proc. R. Soc. Lond. B Biol. Sci.

 

281

:

20140608

.

Google Scholar

Crossref

 

Mitchell

 

J.S.

, 

Rabosky

 

D.L.

 

2016

. 

Bayesian model selection with Bayesian analysis of macroevolutionary mixtures: effects of the model prior on the inferred number of diversification shifts.

 

Methods Ecol. Evol.

 

10.1111/2041-210X.12626

.

Google Scholar

Crossref

 

Moore,

 

B. R.,

 

Höhna,

 

S.,

 

May,

 

M. R.,

 

Rannalla,

 

B.,

 

Huelsenbeck,

 

J. P.

 

2016

. 

Critically evaluating the theory and performance of Bayesian analysis of macroevolutionary mixtures.

 

Proc. Natl. Acad. Sci. USA.

 

113

:

9569

–

9574

.

Google Scholar

Crossref

 

Morlon,

 

H.

 

2014

. 

Phylogenetic approaches for studying diversification.

 

Ecol. Lett.

 

17

:

508

–

525

.

Google Scholar

Crossref

PubMed

 

Morlon

 

H.

, 

Parsons

 

T.L.

, 

Plotkin

 

J.B.

 

2011

. 

Reconciling molecular phylogenies with the fossil record.

 

Proc. Natl. Acad. Sci. USA.

 

108

:

16327

–

16332

.

Google Scholar

Crossref

 

Nee

 

S.

 

2006

. 

Birth-death models in macroevolution.

 

Annu. Rev. Ecol. Evol. Syst.

 

37

:

1

–

17

.

Google Scholar

Crossref

 

Nee

 

S.

, 

May

 

R.M.

, 

Harvey

 

P.H.

 

1994

. 

The reconstructed evolutionary process.

 

Philos. Trans. R. Soc. Lond. B Biol. Sci.

 

344

:

305

–

311

.

Google Scholar

Crossref

PubMed

 

Peters

 

S.E.

 

2005

. 

Geologic constraints on the macroevolutionary history of marine animals.

 

Proc. Natl. Acad. Sci. USA.

 

102

:

12326

–

12331

.

Google Scholar

Crossref

 

Peters

 

S.E.

 

2006

. 

Macrostratografia da América do Norte.

 

Geol.

 

114

:

391-412

.

Google Scholar

Crossref

 

Peters

 

S.E.

, 

Foote

 

M.

 

2001

. 

Biodiversidade no Phanerozic: uma reinterpretação.

 

Paleobiologia.

 

27

:

583-601

.

Google Scholar

Crossref

 

Pujos

 

F.

, 

Iuliis

 

G.D.

, 

Cartelle

 

C.

 

2016

. 

Uma visão geral paleogeográfica das preguiças fósseis tropicais: para uma compreensão da origem das preguiças suspensas?

 

Mamm. O Evol.

 

1-20

.

Google Scholar

 

Purvis

 

A.

 

2008

. 

Abordagens filogenéticas para o estudo da extinção.

 

Annu, annu. Reverendo Ecol. Evol. Syst.

 

39

:

301-319

.

Google Scholar

Crossref

 

Quental

 

T.B.

 

Marshall

 

C.R.

 

2010

. 

Dinâmica da diversidade: filogenias moleculares precisam do registro fóssil.

 

Tendências Ecol. Evol.

 

25

:

434-441

.

Google Scholar

Crossref

PubMed

 

Quental,

 

T.B.,

 

Marshall

 

C. R.

 

2013

. 

Como a Rainha Vermelha leva os mamíferos terrestres à extinção.

 

Ciência.

 

341

:

290-292

.

Google Scholar

Crossref

PubMed

 

Rabosky

 

D.L.

 

2009

. 

A herdabilidade das taxas de extinção associa padrões de diversificação em filogenias moleculares e fósseis.

 

Biol.

 

58

:

629-640

.

Google Scholar

Crossref

PubMed

 

Rabosky

 

D.L.

 

2010

. 

As taxas de extinção não devem ser estimadas a partir de filogenias moleculares.

 

Evolução.

 

64

:

1816-1824

.

Google Scholar

Crossref

PubMed

 

Rabosky

 

D.L.

 

2014

. 

Detecção automática de inovações-chave, mudanças de taxas e dependência da diversidade em árvores filogenéticas.

 

PLos Um.

 

9

:

e89543

.

Google Scholar

Crossref

PubMed

 

Rabosky

 

D.L.

 

2016

. 

Desafios na estimativa de extinção de filogenias moleculares: uma resposta a Beaulieu e O'Meara.

 

Evolução.

 

70

:

218-228

.

Google Scholar

Crossref

PubMed

 

Rabosky

 

D.L.

, 

Grundler

 

M.

, 

Anderson

 

C.

, 

Title

 

P.

, 

Shi

 

J.J.

, 

Brown

 

J.W.

, 

Huang

 

H.

, 

Larson

 

J.G.

 

2014

. 

BAMMtools: um pacote R para a análise da dinâmica evolutiva em árvores filogenéticas.

 

Métodos Ecol. Evol.

 

5

:

701-707

.

Google Scholar

Crossref

 

Rabosky,

 

D.L.,

 

Mitchell,

 

J.S.,

 

Chang

 

J.

 

2017

. 

O BAMM é falho? Preocupações teóricas e práticas na análise de modelos de diversificação multi-taxa.

 

Biol.

 

66

:

477-498

.

Google Scholar

Crossref

PubMed

 

Rabosky

 

D.L.

, 

Santini

 

F.

, 

Eastman

 

J.

, 

Smith

 

S.A.

, 

Sidlauskas

 

B.

, 

Chang

 

J.

, 

Alfaro

 

M.E.

 

2013

. 

As taxas de especiação e evolução morfológica estão correlacionadas entre a maior radiação vertebrada.

 

Nat. Commun.

 

4

:

1958

.

Google Scholar

Crossref

PubMed

 

Raup

 

D.M.

 

Diversidade taxonômica durante o phanerozico.

 

Ciência.

 

177

:

1065-1071

.

Google Scholar

Crossref

PubMed

 

Raup

 

D.M.

 

Diversidade de espécies no Phanerozoico: uma interpretação.

 

Paleobiologia.

 

2

:

289-297

.

Google Scholar

Crossref

 

Raup

 

D.M.

 

Modelos matemáticos de cladogênese.

 

Paleobiologia.

 

11

:

42-52

.

Google Scholar

Crossref

 

Raup

 

D.M.

 

O papel da extinção na evolução.

 

Proc. Natl. Acad. Sci. USA.

 

91

:

6758-6763

.

Google Scholar

Crossref

 

Raup

 

D.M.

, 

Gould

 

S.J.

, 

Schopf

 

T.J.M.

, 

Simberloff

 

D.S.

 

1973

. 

Modelos estocásticos de filogenia e a evolução da diversidade.

 

Geol.

 

81

:

525-542

.

Google Scholar

Crossref

 

Ronquist

 

F.

, 

Klopfstein

 

S.

, 

Vilhelmsen

 

L.

, 

Schulmeister

 

S.

, 

Murray

 

D.L.

, 

Rasnitsyn

 

A.P.

 

2012

. 

A total-evidence approach to dating with fossils, applied to the early radiation of the hymenoptera.

 

Syst. Biol.

 

61

:

973

–

999

.

Google Scholar

Crossref

PubMed

 

Rosenblum

 

E.B.

, 

Sarver

 

B.A.J.

, 

Brown

 

J.W.

, 

Roches

 

S.D.

, 

Hardwick

 

K.M.

, 

Hether

 

T.D.

, 

Eastman

 

J.M.

, 

Pennell

 

M.W.

, 

Harmon

 

L.J.

 

2012

. 

Goldilocks Meets Santa Rosalia: an ephemeral speciation model explains patterns of diversification across time scales.

 

Evol. Biol.

 

39

:

255

–

261

.

Google Scholar

Crossref

PubMed

 

Sakamoto,

 

M.,

 

Benton,

 

M. J.,

 and 

Venditti.

 

C.

 

2016

. 

Dinosaurs in decline tens of millions of years before their final extinction.

 

Proc. Natl. Acad. Sci. USA.

 

113

:

5036

–

5040

.

Google Scholar

Crossref

 

Sepkoski

 

J.J.

 

1978

. 

A kinetic model of phanerozoic taxonomic diversity I.

 

Analysis of marine orders. Paleobiology.

 

4

:

223

–

251

.

Google Scholar

Crossref

 

Sepkoski

 

J.J.

, 

Bambach

 

R.K.

, 

Raup

 

D.M.

, 

Valentine

 

J.W.

 

1981

. 

Phanerozoic marine diversity and the fossil record.

 

Nature.

 

293

:

435

–

437

.

Google Scholar

Crossref

 

Sheets

 

H.D.

, 

Mitchell

 

C.E.

, 

Melchin

 

M.J.

, 

Loxton

 

J.

, 

Štorch

 

P.

, 

Carlucci

 

K.L.

, 

Hawkins

 

A.D.

 

2016

. 

Graptolite community responses to global climate change and the Late Ordovician mass extinction.

 

Proc. Natl. Acad. Sci. USA.

 

113

:

8380

–

8385

.

Google Scholar

Crossref

 

Silvestro

 

D.

, 

Schnitzler

 

J.

, 

Liow

 

L.H.

, 

Antonelli

 

A.

, 

Salamin

 

N.

 

2014a

. 

Bayesian estimation of speciation and extinction from incomplete fossil occurrence data.

 

Syst. Biol.

 

63

:

349

–

367

.

Google Scholar

Crossref

 

Silvestro

 

D.

, 

Salamin

 

N.

, 

Schnitzler

 

J.

 

2014b

. 

PyRate: a new program to estimate speciation and extinction rates from incomplete fossil data.

 

Methods Ecol. Evol.

 

5

:

1126

–

1131

.

Google Scholar

Crossref

 

Stadler

 

T.

 

2010

. 

Sampling-through-time in birth–death trees.

 

J. Theor. Biol.

 

267

:

396

–

404

.

Google Scholar

Crossref

PubMed

 

Stadler

 

T.

 

2013

. 

Recovering speciation and extinction dynamics based on phylogenies.

 

J. Evol. Biol.

 

26

:

1203

–

1219

.

Google Scholar

Crossref

PubMed

 

Stanley

 

S.M.

 

Macroevolução, padrão e processo

. 

Baltimore, Maryland

: 

Johns Hopkins University Press

.

Google Scholar

Stock

 

C.

 

1929

. 

Um censo dos mamíferos pleistocenos do Rancho La Brea, baseado nas coleções do Museu de Los Angeles.

 

Mamífero.

 

10

:

281-289

.

Google Scholar

Crossref

 

Valentine

 

J.W.

 

1989

. 

Quão bom foi o registro fóssil? Pistas do Pleistoceno californiano.

 

Paleobiologia.

 

15

:

83-94

.

Google Scholar

Crossref

 

Valentine

 

J.W.

, 

Jablonski

 

D.

, 

Kidwell

 

S.

, 

Roy

 

K.

 

2006

. 

Avaliando a fidelidade do registro fóssil usando biválvulas marinhas.

 

Proc. Natl. Acad. Sci. USA.

 

103

:

6599-6604

.

Google Scholar

Crossref

 

Varela

 

S.

, 

González-Hernández

 

J.

, 

Sgarbi

 

L.F.

, 

Marshall

 

C.

, 

Uhen

 

M.D.

, 

Peters

 

S.

, 

McClennen

 

M.

 

2015

. 

paleobioDB: um pacote R para baixar, visualizar e processar dados do Banco de Dados de Paleobiologia.

 

Ecografia.

 

38

:

419-425

.

Google Scholar

Crossref

 

Zanno

 

L.E.

, 

Makovicky

 

P.J.

 

2011

. 

Ecomorfologia herbívora e padrões de especialização na evolução dos dinossauros terópodes.

 

Proc. Natl. Acad. Sci. USA.

 

108

:

232-237

.

Google Scholar

Crossref

 

Zhang

 

C.

, 

Stadler

 

T.

, 

Klopfstein

 

S.

, 

Heath

 

T.A.

, 

Ronquist

 

F.

 

2016

. 

Evidências totais datando sob o processo fossilizado de nascimento e morte.

 

Biol.

 

65

:

228-249

.

Google Scholar

Crossref

PubMed

 

© O Autor(s) 2018. Publicado pela Oxford University Press, em nome da Sociedade de Biólogos Sistemáticos. Todos os direitos reservados. Para permissões, envie um e-mail: journals.permissions@oup.com

Este artigo é publicado e distribuído sob os termos da Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)