• Diversificación: tasa a la cual se acumulan los linajes en una filogenia
• Incluye pero no se limita a cladogenesis
• Diversificación: especiación + extinción
• Se ha aplicado en la evolución de especies, poblaciones y lenguajes

Incluye:

• Tasas de especiación y extinción
• Diversificación diferencial entre clados
• Cambios en diversificación asociado a regímenes selectivos
• Covariación entre tasas de diversificación y tasa evolutiva de rasgos

Modelos generalmente incluyen 2 parámetros:

• b = Especiación (birth)
• d = Extinción (death)

Modelo de Yule:

• Un único parámetro (b)

library(phytools)

#simular arbol solo con b 
set.seed(20)
arb1 <- pbtree(n = 100)
plotTree(arb1, ftype = "off")

## Loading required package: ape

## Loading required package: maps

## 
##  # maps v3.1: updated 'world': all lakes moved to separate new #
##  # 'lakes' database. Type '?world' or 'news(package="maps")'.  #

Modelo de Yule:

• Un único parámetro (b)

Modelo especiación-extinción (birth-death):

• 2 parámetros (b & d)

#simular arbol con b y d
set.seed(46)
arb2 <- pbtree(n = 100, b = 0.9, d =  0.4, scale = T)
plotTree(arb2, ftype = "off")

Modelo especiación-extinción (birth-death):

Modelo especiación-extinción (birth-death):

• Acumulación de especies en el tiempo

#simular arbol con b y d
par(mar= c(12,4,0,4))
lttbd <- ltt(arb2)

¿Como inferir tasas de especiación/extinción?

• La extinción deja una "firma" detectable en las filogenias

#remover las especies extintas
par(mfrow = c(1,2))
plotTree(arb2, ftype = "off", mar= rep(1, 4))
title("Con spp extintas")

arb2.2 <- drop.tip(arb2, getExtinct(arb2))
plotTree(arb2.2, ftype = "off", mar= rep(1, 4))
title("Sin spp extintas")

¿Como inferir tasas de especiación/extinción?

• La extinción deja una "firma" detectable en las filogenias

¿Como inferir tasas de especiación/extinción?

• La extinción deja una "firma" detectable en las filogenias

par(mfrow = c(1,1))

#curvas de acumulacion de spp
lttbd <- ltt(arb2.2)
lttbd2 <- ltt(arb2, plot = F)
lines(lttbd2$times, log(lttbd2$ltt), col = "red")

¿Como inferir tasas de especiación/extinción?

• La extinción deja una "firma" detectable en las filogenias

¿Como inferir tasas de especiación/extinción?

• La extinción deja una "firma" detectable en las filogenias

#curvas de acumulacion de spp
lttbd <- ltt(arb2.2)
lines(c(0, max(nodeHeights(arb2.2))), c(log(2), 
      log(length(arb2.2$tip.label))), 
    lty = "dashed", lwd = 2, col = "red")

¿Como inferir tasas de especiación/extinción?

• La extinción deja una "firma" detectable en las filogenias

¿Como inferir tasas de especiación/extinción?

• La extinción deja una "firma" detectable en las filogenias

#arebol coalescente
arbcoal <- rcoal(n = 100)
plotTree(arbcoal, ftype = "off", mar= rep(1, 4))

¿Como inferir tasas de especiación/extinción?

• La extinción deja una "firma" detectable en las filogenias

lttac <- ltt(arbcoal)
lines(c(0, max(nodeHeights(arbcoal))), c(log(2), 
      log(length(arbcoal$tip.label))), 
    lty = "dashed", lwd = 2, col = "red")

• b - d = tasa neta de diversificación (acumulación nuevas especies)
• d / b = tasa de extinción relativa (~la vulnerabilidad de extinción de clados pequeños)

#install.packages("laser")
library(laser)

## Loading required package: geiger

tiemp.ram <- branching.times(arb2.2)
        param.bd <- bd(tiemp.ram) 

Tasa neta de diversificación del árbol

0.5/0.9

## [1] 0.5555556

Tasa neta de diversificación estimada

param.bd$r1

## [1] 3.809059

Tasa de extinción relativa del árbol

0.9 - 0.5

## [1] 0.4

Tasa de extinción relativa estimada

param.bd$a

## [1] 0.4198187

Incertidumbre

• Estimaciones son mas precisas cuando se incluyen todas las especies (incluyendo extintas)
• d es generalmente mas impreciso que b
• Estimaciones son poco confiables cuando se violan supuestos (e.g. tasas constante de diversificación)
• Son mas confiables en términos relativos (e.g. comparación entre clados)

Algunos supuestos pueden ser puestos a prueba:

• Tasa constante de diversificación

gammatest(lttbd)

## $gamma
## [1] 1.425982
## 
## $p
## [1] 0.1538736

Ejemplo con filogenia de Psitacidos (Arini)

arini.datos.comp <-
readRDS(gzcon(url("http://marceloarayasalas.weebly.com/uploads/2/5/5/2/25524573/arini.datos.comp.rds")))

#definir arbol y datos
arb.arini <- arini.datos.comp[[1]]

plotTree(arb.arini, fsize = 0.6)

Ejemplo con filogenia de Psitacidos (Arini)

tiemp.ram <- branching.times(arb.arini)

bd(tiemp.ram) 

## $LH
##     LH.52 
## -23.57133 
## 
## $r1
##    smax.52 
## 0.08117029 
## 
## $a
## [1] 0
## 
## $aic
##    LH.52 
## 51.14265

Ejemplo con filogenia de Psitacidos (Arini)

lttArini <- ltt(arb.arini, plot = F)
gammatest(lttArini)

## $gamma
## [1] -0.8740852
## 
## $p
## [1] 0.3820719

Diversificación varia en el tiempo o entre clados:

• Cambios en diversificación pueden estar asociados a la evolución de nuevos estados en un rasgo
• Cambios en diversificación pueden estar asociados a la variación de un rasgo continuo

Cambios asociados a estados en un rasgo binario

• Modelo BiSSE (Binary state speciation and extinction, Madisson et al. 2007)

Fuente: Goldberg, E. E., Lancaster, L. T., & Ree, R. H. (2011). Phylogenetic inference of reciprocal effects between geographic range evolution and diversification. Systematic Biology, 60(4), 451-465.

Simulación de evolucion bajo el Modelo BiSSE:

• Simular un rasgo binario sobre una filogenia
• tasas de transicion entre estados
• c(λ0, λ1, μ0, μ1, q01, q10) = c(0.1, 0.3, 0.03, 0.03, 0.01, 0.01)

Simulación de evolucion bajo el Modelo BiSSE:

library(diversitree)

set.seed(75)
tasas <- c(0.1, 0.3, 0.03, 0.03, 0.01, 0.01)
arb.bisse <- tree.bisse(pars = tasas, max.t=30, x0=0)
colors <- c("maroon", "blue")
plot(history.from.sim.discrete(arb.bisse, 0:1), arb.bisse,
     col=colors, show.tip.label=FALSE)

Simulación de evolucion bajo el Modelo BiSSE:

Ajuste del modelo BiSSE:

• Estima si la diversificacion es dependiente del estado del rasgo binario
• Compara el ajuste de un modelo donde las tasas de transicion pueden variar contra uno donde son restringidos

Ajuste del modelo BiSSE: - Modelo con tasas no restringidas tasas

datos.bisse <- make.bisse(arb.bisse, arb.bisse$tip.state)

p <- starting.point.bisse(arb.bisse)
ajuste.no.restr <- find.mle(datos.bisse, p)

coef(ajuste.no.restr)

##      lambda0      lambda1          mu0          mu1          q01 
## 5.810208e-02 3.782662e-01 9.558864e-07 1.192523e-01 5.033573e-08 
##          q10 
## 1.747444e-02

Ajuste del modelo BiSSE:

# modelo donde varian las tasas
bisse.restr = constrain(datos.bisse, lambda0 ~ lambda1, mu0 ~ mu1)

ajuste.restr = find.mle(bisse.restr, p)

## Warning in guess.constrained.start(func, x.init): Guessing parameters while
## constraining model - may do badly

coef(ajuste.restr)

##      lambda1          mu1          q01          q10 
## 3.432876e-01 1.361948e-01 1.818086e-08 2.115162e-02

Ajuste del modelo BiSSE:

AIC(ajuste.no.restr, ajuste.restr)

##                 df      AIC
## ajuste.no.restr  6 772.1837
## ajuste.restr     4 792.6644

anova(ajuste.no.restr, ajuste.restr)

##         Df   lnLik    AIC  ChiSq Pr(>|Chi|)    
## full     6 -380.09 772.18                      
## model 1  4 -392.33 792.66 24.481  4.832e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ajuste del modelo BiSSE sobre datos Arini:

datos.arini <- arini.datos.comp[[2]]
arb.arini <- chronoMPL(arb.arini)
is.ultrametric(arb.arini)

## [1] TRUE

veg.bin <- cut((datos.arini$Vegetation.index.EVI), breaks = 2)

#renombrar estados
levels(veg.bin) <- c(0, 1)
veg.bin <- as.numeric(as.character(veg.bin))

names(veg.bin) <- datos.arini$Species

Ajuste del modelo BiSSE sobre datos Arini:

#generar funcion LL para BiSSE
arini.bisse <- make.bisse(arb.arini, veg.bin)

#correr modelo
p <- starting.point.bisse(arb.arini)
arini.no.restr <- find.mle(arini.bisse, p)

#Parametros estimados
coef(arini.no.restr)

##      lambda0      lambda1          mu0          mu1          q01 
## 1.667172e-01 4.880700e-02 9.518747e-06 5.037510e-06 3.540355e-01 
##          q10 
## 9.398720e-02

Ajuste del modelo BiSSE sobre datos Arini:

#generar funcion LL para BiSSE
arini.restr = constrain(arini.bisse, lambda0 ~ lambda1, mu0 ~ mu1)

#correr modelo
arini.restr = find.mle(arini.restr, p)

## Warning in guess.constrained.start(func, x.init): Guessing parameters while
## constraining model - may do badly

coef(arini.restr)

##      lambda1          mu1          q01          q10 
## 8.118031e-02 2.264931e-07 1.964068e-01 5.630865e-02

Ajuste del modelo BiSSE sobre datos Arini:

#comparar modelos
anova(arini.restr, arini.no.restr)

##         Df   lnLik    AIC ChiSq Pr(>|Chi|)
## minimal  4 -198.65 405.30                 
## model 1  6 -197.78 407.56 1.744     0.4181

Modelo MuSSE (Multi-state speciation and extinction, Pagel 1994)

• Cambios asociados a estados en un rasgo multi-estado
• Similar a BiSSE
• Mas parametros estimados -> computacionalmente intensivo

Simulación de evolucion bajo el Modelo BiSSE:

• Simular un rasgo multi-estado sobre una filogenia
• Solo transiciones entre estados adyacentes (1 <-> 2 <-> 3)
• Especiacion y extincion aumentan con el valor del estado (1 -> 2 -> 3)
• Tasas de transicion entre estados:
  • c(λ1, λ2, λ3, μ1, μ2, μ3, q12, q13, q21, q23, q31, q32) = c(0.1, 0.15, 0.2, 0.03, 0.045, 0.06, 0.05, 0, 0.05, 0.05, 0, 0.05)

Simulación de evolucion bajo el Modelo BiSSE:

set.seed(2)
datos.muse <- tree.musse(pars = c(0.1, 0.15, 0.2, 0.03, 0.045, 0.06, 0.05,
                                  0, 0.05, 0.05, 0, 0.05), 30, x0=1)

# graficar rasgo + filogenia
col <- c("blue", "orange", "red")
h <- history.from.sim.discrete(datos.muse, 1:3)
plot(h, datos.muse, cex=.7, col=col, show.tip.label = F)

Simulación de evolucion bajo el Modelo BiSSE:

Simulación de evolucion bajo el Modelo BiSSE:

estados <- datos.muse$tip.state

#funcion LL
lik <- make.musse(datos.muse, estados, 3)

#modelo restringido
musse.restr <- constrain(lik, lambda2 ~ lambda1, lambda3 ~ lambda1,
mu2 ~ mu1, mu3 ~ mu1, q13 ~ 0, q21 ~ q12, q23 ~ q12, q31 ~ 0, q32 ~ q12)

#busqueda heurística
p <- starting.point.musse(datos.muse, 3)
ajuste.restr <- find.mle(musse.restr, p[argnames(musse.restr)])

Simulación de evolucion bajo el Modelo BiSSE:

# Modelo donde especiación varia
musse.norestr <- constrain(lik, mu2 ~ mu1, mu3 ~ mu1,
q13 ~ 0, q21 ~ q12, q23 ~ q12, q31 ~ 0, q32 ~ q12)

ajuste.norestr <- find.mle(musse.norestr, p[argnames(musse.norestr)])

anova(ajuste.norestr, ajuste.restr)

##         Df   lnLik    AIC   ChiSq Pr(>|Chi|)
## full     5 -110.21 230.42                   
## model 1  3 -110.64 227.27 0.84603     0.6551

1. Simule un proceso evolutivo con divergencia dependiente de un rasgo binario (BiSSE) donde las tasas de especiación para los estados del rasgo sean las mismas pero las tasas de extinción sean 4 veces mayores en un estado.

2. Compare el ajuste de un modelo no-restringido y uno restringido sobre el proceso simulado.

3. Simule un proceso evolutivo con divergencia dependiente de un rasgo con 3 estados (MuSSE) sobre los datos de Arini.