Patrones taxonómicos, geográficos y filogenéticos
Parte I
Juan D. Vásquez-Restrepo
Department of Ecology & Evolutionary Biology, Princeton Universitydvasquez@princeton.edu
Cuando realizamos análisis a grandes escalas, como los que son habituales en macroecología y macroevolución, es importante reconocer que el concepto de escala puede adoptar distintas formas. No solo implica trabajar con un gran número de especies o con grupos taxonómicos amplios, sino también abordar regiones geográficas extensas que abarcan múltiples contextos ecológicos.
Recientemente publicamos un artículo en el que describimos los patrones taxonómicos, geográficos y filogenéticos en el estado de conservación de los reptiles escamados de Colombia (Vásquez-Restrepo & García-Cobos 2026). Si bien los productos se incluyeron como material suplementario, considero importante mostrar cómo se calculan las distintas métricas utilizadas y cómo estas pueden representarse espacialmente.
El objetivo principal de este tutorial no es profundizar en la interpretación de los índices o resultados, sino centrarse en el procesamiento de los datos. Este paso suele ser el más desafiante para quienes buscan reproducir este tipo de análisis o adaptarlos a otros sistemas de estudio. Al igual que en otros tutoriales que he publicado, la intención es presentar los contenidos de forma clara, progresiva y didáctica, para que sean accesibles a personas con distintos niveles de experiencia en programación.
1. Antes de comenzar
Para este tutorial, dado que exploraremos patrones taxonómicos, geográficos y filogenéticos, es necesario contar con tres insumos básicos: una lista de especies, información espacial asociada a esas especies y un árbol que represente sus relaciones evolutivas. Los datos que vamos a utilizar pueden descargarse directamente desde mi GitHub.
Para facilitar el flujo de trabajo, los datos han sido preprocesados de modo que sabemos cómo aparece cada especie en los distintos conjuntos de datos. Para ello, contamos con un par de columnas adicionales que funcionarán como un diccionario taxonómico. Este punto es particularmente importante, ya que muchos análisis requieren hacer coincidir estos objetos entre sí y los nombres de las especies pueden variar entre fuentes por distintas razones, como cambios taxonómicos o errores tipográficos. Estas inconsistencias pueden generar errores en los análisis o provocar la exclusión involuntaria de especies. Si bien no es estrictamente necesario que el 100% de las especies incluidas en el árbol tengan información geográfica (o viceversa), es deseable que una proporción alta esté representada en ambos conjuntos de datos. De este modo, se reduce el sesgo de muestreo y se mejora la robustez de los patrones que se buscan identificar.
En cuanto a los datos de distribución, es posible utilizar los polígonos de la UICN, que para este ejercicio corresponden a la versión 2023. Mientras que, para las relaciones filogenéticas, existen varias alternativas: generar un árbol propio, utilizar árboles publicados, o descargar filogenias disponibles en repositorios públicos como VertLife. En términos prácticos, para la mayoría de los análisis basta con un árbol filogenético dicotómico. No obstante, también es posible trabajar con un conjunto de árboles para incorporar la incertidumbre filogenética (Rangel et al. 2015). En este ejemplo utilizaremos una muestra de 100 árboles, ya que no contamos con un árbol completamente bifurcado que incluya a la mayoría de nuestras especies de interés. Estos árboles provienen de una distribución de probabilidades (Tonini et al. 2015) y permiten capturar parte de esa incertidumbre en los análisis posteriores.
Importando los datos
# Cargamos las librerías necesarias
# ---Manejo de datos y visualización general
library(tidyverse)
library(cowplot)
library(biscale)
library(ggtree) # Nota: Instalar con remotes::install_github('YuLab-SMU/ggtree')
# ---Análisis espacial (GIS)
library(sf)
library(terra)
library(raster)
library(nngeo)
library(wdpar)
# --- Filogenética
library(ape)
library(phytools)
library(caper)
# ---Macroecología y distribución
library(letsR)
library(speciesRaster) # Nota: Instalar con remotes::install_github('ptitle/speciesRaster')
# ---Filogenética espacial y diversidad
library(picante)
library(phyloregion)
library(PhyloMeasures) # Nota: Instalar con remotes::install_github('cran/PhyloMeasures')
# ---Utilidades matemáticas
library(matrixStats)
# Definimos nuestra ruta de trabajo
setwd("D:/Daniel/Desktop/resources")
# Generamos una paleta de color para los mapas (nos será útil después)
pal <- rev(terrain.colors(100))
# Cargamos nuestra lista de especies
species <- read.csv("species_list.csv")
# Cargamos nuestro conjunto de árboles
trees <- read.tree("trees/Tonini_trees.trees")
# Cargamos los polígonos de distribución
IUCN <- st_read(paste0(getwd(), "/geo"), "IUCN_2023.1_Colombia", quiet = TRUE)Cuando trabajamos con datos geográficos, es fundamental verificar el sistema de coordenadas, ya que un sistema faltante, incorrecto o inconsistente puede convertirse en un problema serio en etapas posteriores. Existen distintos sistemas de referencia espacial y su elección depende del objetivo del análisis. De forma general, pueden agruparse en sistemas planos (o proyectados) y sistemas geográficos.
Aquí no profundizaré en las diferencias entre ellos ni en los criterios para elegir uno u otro, ya que es un tema que merece un tratamiento más detallado y queda fuera del alcance de este tutorial (tarea para casa). Para nuestros análisis utilizaremos coordenadas geográficas expresadas como longitud y latitud, a partir de un estándar ampliamente conocido: WGS84. Los datos provistos ya se encuentran en este sistema, lo cual podemos verificar fácilmente en R.
# Validamos el sistema de referencia (debe indicar WGS84)
st_crs(IUCN)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]# Si el sistema de coordenadas es "unknown" o "NULL" pero sabemos que debe ser WGS84, podemos asignarlo manualmente
st_crs(IUCN) <- 4326 # Correr SOLO si es el caso
# Si por el contrario el sistema de coordenadas es diferente, debemos transformarlo
IUCN <- st_transform(IUCN, crs = 4326) # Correr SOLO si es el casoValidaciones iniciales
Una vez que tenemos los datos cargados, es necesario aplicar algunos filtros y transformaciones antes de continuar. Para este tutorial, los polígonos de distribución no solo han sido recortados, sino también filtrados y corregidos previamente, con el objetivo de reducir el tamaño de los archivos y facilitar su manejo. En contraste, los árboles filogenéticos no han sido modificados. Por esta razón, es necesario reducir el conjunto de especies presentes en los árboles para que coincida con aquellas incluidas en nuestra lista base.
Aquí es donde haremos uso de nuestro diccionario taxonómico. Personalmente no me gusta modificar el nombre de los polígonos o los terminales del árbol directamente, por lo que prefiero agregar a mi tabla base un par de columnas que indiquen el nombre con el que aparece cada especie en cada conjunto de datos. En este caso, utilizaremos “Polygon.name” y “Phylo.name” para los polígonos y el árbol, respectivamente. Realizar este cruce puede ser una tarea tediosa, especialmente cuando los datos provienen de años diferentes. Como se mencionó antes, la taxonomía es cambiante y cuanto mayor sea la distancia temporal entre los conjuntos de datos, mayor será la probabilidad de tener que lidiar con sinónimos o cambios taxonómicos.
# Validamos la clase de nuestros polígonos (debe ser "sf" y "data.frame")
class(IUCN)## [1] "sf" "data.frame"# Convertimos los polígonos a un objeto espacial de R
IUCN <- as(IUCN, "Spatial")
# Validamos de nuevo la clase de los polígonos (ahora debe ser "SpatialPolygonsDataFrame")
class(IUCN)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"# Filtramos nuestra lista de especies para eliminar aquellas que no están en filogenia
tips <- species$Phylo.name[!is.na(species$Phylo.name)]
# Validamos cuántas especies tenemos mapeadas en filogenia (debe ser mayor a cero)
length(tips)## [1] 551# Filtramos nuestros árboles para incluir solo las especies de nuestra lista filtrada
trees_pruned <- keep.tip(trees, tips)
# Validamos cuántas especies hay en los árboles (debe dar igual al número mapeado en filogenia)
# ---Como tenemos una lísta y no un único árbol, indexamos cualquiera (con verificar el primero es suficiente)
length(trees_pruned[[1]]$tip.label)## [1] 551En el caso de los polígonos de distribución, es importante aclarar que, en ocasiones, una misma especie puede estar representada por múltiples features (o polígonos). Por ejemplo, la UICN suele incluir varios polígonos para una misma especie, dependiendo de cómo se categoricen sus áreas de distribución, por ejemplo, si son áreas de residencia o de introducción.
Para poder realizar los análisis de este tutorial, es necesario
contar con un único polígono por especie. Esto implica eliminar aquellas
áreas que no sean de interés para el análisis o alternativamente,
fusionar los distintos polígonos en un solo feature. Aunque
este procedimiento puede realizarse directamente en R, suele ser más
simple y rápido usando QGIS (cualquier software GIS sirve, pero este es
Open Source), por ejemplo, mediante la función Disolver.
Adicionalmente, algunos polígonos pueden presentar geometrías inválidas,
lo que puede generar errores en etapas posteriores. Por esta razón,
también es necesario corregir las geometrías. Esto último puede hacerse
en QGIS mediante la herramienta Reparar geometrías.
Si las geometrías son inválidas, nos daremos cuenta en el siguiente paso, ya que no será posible continuar con el análisis. En cambio, si existe más de un feature por especie, esto no será evidente a menos que se haya validado previamente. Por esta razón, vamos a comparar el número de polígonos y el número de nombres únicos asociados a ellos. De no hacerlo, podríamos terminar con especies duplicadas que, para el análisis, se contabilizan como entidades distintas.
# Validamos cuántos polígonos tenemos
length(IUCN)## [1] 486# Validamos cuántas especies únicas hay en los polígonos (debe dar igual an anterior)
length(unique(IUCN@data$BINOMIAL))## [1] 4862. Diversidad α espacial
Ahora que tenemos nuestros datos filtrados y curados, es hora de poner manos a la obra. Lo primero será calcular la diversidad α espacial. En un tutorial anterior ya expliqué qué es la diversidad α, que en este contexto no es más que el número de especies dentro de cada comunidad. Ahora bien, ¿cuáles son nuestras comunidades?
Este tipo de análisis suele realizarse superponiendo una cuadrícula sobre el área de interés. En otras palabras, dividimos el espacio en celdas de un tamaño determinado y cada una de ellas se considera una comunidad. Para el cálculo de la diversidad α, lo que hacemos es contar el número de polígonos que se sobrelapan en cada celda. Por simplicidad, trabajaremos con celdas cuadradas, aunque no es la única opción, también es posible utilizar hexágonos (de hecho son más adecuados). La elección entre una u otra forma dependerá del sistema de referencia espacial que estemos usando y del objetivo del análisis.
Es importante tener presente que la forma en que representamos los mapas en un plano está sesgada, ya que las áreas en latitudes altas suelen verse distorsionadas. Por ejemplo, en el ecuador un grado equivale a aproximadamente 111 km, mientras que en la latitud 40° esta distancia se reduce a unos 85 km. Como estamos trabajando a una escala espacial relativamente grande, podemos definir nuestro grano entre uno y medio grado. Es común querer realizar análisis muy detallados, pero siempre es necesario considerar el costo-beneficio. Por ejemplo, aunque es posible obtener datos espaciales a una resolución de 30 metros, estos pueden ser tan pesados y el análisis tan costoso computacionalmente que no solo podría tomar días (o meses), sino que los patrones observados podrían no diferir de aquellos obtenidos con una resolución más gruesa. Aquí es donde se vuelve importante entender nuestros datos. Los polígonos de la UICN que estamos utilizando son gruesos, por lo que los errores de comisión, es decir, la sobrepredicción, tienden a aumentar a medida que reducimos la escala.
# Definimos la extensión del análisis
# ---En este caso, las coordenadas corresponden a una caja que abarca todo Colombia
ex <- extent(c(-84, -64, -6, 14))
# Validamos la existencia de una columna llamada "BINOMIAL" en nuestros datos espaciales
# ---El paquete "letsR" es quisquilloso, si no existe esta columna, no corre
"BINOMIAL" %in% colnames(IUCN@data) # Debe devolver TRUE## [1] TRUE# Calculamos la matriz de presencia-ausencia a partir de los polígonos
# ---Aquí es donde necesitamos la columna llamada "BINOMIAL"
paMa <- lets.presab(IUCN,
xmn = bbox(ex)[1,1], # X mínimo
xmx = bbox(ex)[1,2], # X máximo
ymn = bbox(ex)[2,1], # Y mínimo
ymx = bbox(ex)[2,2], # Y Máximo
resol = 0.5, # Resolución espacial en unidades del mapa (0.5°)
crs = crs(IUCN), # Sistema de referencia general (hereda WGS84)
crs.grid = crs(IUCN), # Sistema de referencia de la grilla (hereda WGS84)
count = FALSE) # Mostrar el progreso (aquí uso FALSE para no saturar la consola, pero idealmente TRUE)Con esto quedaría listo el primer paso. Dentro del objeto
paMa tenemos una matriz de presencia-ausencia y su
correspondiente raster, a los cuales podemos acceder fácilmente. Aunque
para análisis futuros podríamos simplemente reutilizar el objeto que
acabamos de crear, por experiencia propia he notado que, incluso si lo
exportamos como un archivo .RData, al volver a abrir la
sesión suele generar errores y obliga a recalcularlo. Este problema está
relacionado con la transición del paquete raster a
terra (al mejorar unas cosas, arruinaron otras).
Durante este cambio, los desarrolladores modificaron la forma como se
procesan los rasters espaciales, pasando de una clase
Raster a SpatRast. Esto tiene varias ventajas,
ya que es mucho más rápido en análisis con grandes volúmenes de píxeles
pues no carga todo en la memoria RAM. Sin embargo, ahora los raster no
existen en la sesión sino en un archivo temporal enlazado, por lo que al
cargar objetos generados en sesiones anteriores, estos simplemente se
rompen. Por esta razón, tenemos que convertir nuestro raster en un
objeto autosuficiente antes de continuar.
# Volvemos el raster un objeto autosuficiente para que no se rompa
# ---NO hacer con rasters pesados o gigantes (en ese caso, es mejor exportarlos como .tiff)
paMa$Richness_Raster <- terra::wrap(paMa$Richness_Raster)
# Graficamos el raster
plot(terra::unwrap(paMa$Richness_Raster), col = pal)
3. Diversidad β espacial
Vamos a continuar ahora con la diversidad β. Desafortunadamente, aún no tengo un tutorial en línea que explique algunos de los fundamentos teóricos de esta dimensión de la biodiversidad. Sin embargo, cuento con un borrador de capítulo de libro, todavía sin revisar, que pongo a disposición del lector. En este contexto, entenderemos la diversidad β espacial como recambio (aunque también es posible calcular el anidamiento). El recambio no es más que una forma de cuantificar qué tan parecidas o diferentes son dos comunidades. En la ecología tradicional, solemos comparar una comunidad A con una comunidad B. En ecología espacial, en cambio, nuestras comunidades pueden ser cientos o miles, tantas como píxeles o celdas en nuestra matriz.
Por esta razón, es necesario abordar el problema con un enfoque ligeramente diferente, ya que no tiene sentido comparar todas las comunidades entre sí. Por ejemplo, el recambio entre un píxel ubicado en la Alta Guajira (un desierto) y otro en el extremo sur de la Amazonía colombiana (selva tropical), separados por más de 1700 km, no solo será total, sino que además no aporta información útil. En este sentido, lo que haremos será comparar comunidades cercanas entre sí. Para ello, introduciremos el concepto de ventana móvil. Una ventana móvil es, literalmente, una ventana que se desplaza a través de nuestra matriz, generando pequeños subgrupos sobre los cuales realizaremos los cálculos. En otras palabras, para cada píxel se define un vecindario, se calcula la diversidad β y luego se pasa al siguiente píxel, donde el proceso se repite.
Lo primero que tenemos que hacer es generar una serie de rasters apilados (o stack), uno por cada especie. Esto lo podemos lograr fácilmente a partir de la matriz de presencia-ausencia que generamos en el paso anterior. Esta matriz contiene la información de longitud y latitud de cada celda en las dos primeras columnas y la presencia-ausencia de cada especie a partir de la tercera. De este modo, podemos recorrer la matriz con un bucle, extrayendo la información necesaria para cada especie. Este procedimiento es necesario porque el paquete que utilizaremos requiere los datos en este formato.
# Generamos un stack vacío
stack <- stack()
# Recorremos la matriz de presencia-ausencia para rellenar nuestro stack
# ---La lógica es simple, lon/lat es igual para todos solo cambia la columna de la especie
for(i in 3:ncol(paMa$Presence_and_Absence_Matrix)){
tmp.ras <- rasterFromXYZ(paMa$Presence_and_Absence_Matrix[, c(1, 2, i)])
stack <- addLayer(stack, tmp.ras)
}
# Convertimos el stack a un objeto de tipo "SpeciesRaster"
speRas <- createSpeciesRaster(stack)
# Calculamos la diversidad beta espacial
beta <- betaDiversity_speciesRaster(speRas, radius = 3, metric = "sorensen")
# Asignamos cero a los NA
beta[is.na(beta)] <- 0Para el procedimiento anterior utilizamos el paquete
speciesRaster, que si revisamos su página de GitHub, veremos que
se presenta como el precursor de otro paquete más reciente llamado
epm. La razón por la cual usamos el paquete antiguo es que,
en mi experiencia, resulta más intuitivo. Por ejemplo, el paquete
speciesRaster recibe como argumento el radio de la ventana
expresado en número de celdas. Así, si queremos una ventana de 25 celdas
cuadradas (5 × 5), necesitamos un radio de dos, es decir, dos celdas en
cada dirección desde el centro. Por el contrario, epm
recibe el radio de la ventana en unidades del mapa, lo que implicaría
convertir el número de celdas a grados.
Por su parte, epm permite calcular la diversidad β
particionada (Baselga 2010), lo cual es super útil en ciertos contextos
donde nos interesa el recambio + anidamiento. Sin embargo, aquí estamos
utilizando el recambio combinado ya que tiene una interpretación más
directa. Así, una comunidad con un valor de 0.8 tendrá una diferencia de
80% en su composición respecto a otra. En cualquier caso, ambos paquetes
deberían dar lo mismo. El problema es que, si en epm nos
dejamos llevar por el nombre de los componentes y calculamos el solo
recambio, los valores obtenidos serán distintos a los que ya calculamos
(unas diez veces más pequeños). Esto ocurre porque, al descomponer la
diversidad β, el recambio queda libre de la influencia de las
diferencias en riqueza. Para obtener resultados comparables, sería
necesario calcular los componentes aditivos, es decir, recambio +
anidamiento.
Una vez hecho esto, solo nos queda ver el resultado. Para mantener un
estilo visual consistente, a partir de ahora transformaremos todos los
rasters en objetos SpatRast, ya que su estilo visual por
defecto es más agradable que el estilo base de R.
# Convertimos el raster de diversidad beta a un objeto SpatRas autosuficiente
beta <- as(beta, "SpatRaster")
beta <- terra::wrap(beta)
# Graficamos el raster
plot(terra::unwrap(beta), col = pal)
Antes de continuar quisiera mencionar que, en el raster anterior es posible observar valores muy altos de recambio en la región del Pacífico. Esto ocurre porque allí se registra una sola especie de escamado marino (Hydrophis platurus). Si conocemos nuestros datos y entendemos cómo funciona la diversidad β, podemos deducir que este patrón es un artificio. No tiene mucho sentido comparar el recambio con comunidades que contienen una sola especie que, además, al ser marina, no se encuentra en comunidades terrestres.
Por esta razón, en nuestro artículo mostramos este caso de manera separada. Esto puede lograrse fácilmente recortando el área terrestre de la marina y graficando cada ecosistema por separado. Otra opción es superponer ambos resultados, pero es importante tener en cuenta que los valores máximos de cada región serán distintos, por lo que esto debe aclararse explícitamente o representarse usando dos escalas diferentes.
# Cargamos el contorno del América
Ame <- st_read(paste0(getwd(), "/geo"), "America_MERGED", quiet = TRUE)
# Enmascaramos la zona terrestre para la diversidad beta
beta_land <- mask(terra::unwrap(beta), Ame)
# Enmascaramos la zona marina para la diversidad beta
beta_sea <- mask(terra::unwrap(beta), Ame, inverse = TRUE)
# Graficamos cada segmento
plot(beta_land, col = pal)
plot(beta_sea, col = pal)
4. Especies endémicas y amenazadas
Otro de los patrones espaciales que analizamos en nuestro trabajo fue la distribución geográfica de las especies endémicas y amenazadas. Para ello, solo necesitamos filtrar las especies de interés tanto en nuestra tabla base como en la matriz de presencia-ausencia. Este paso es sencillo, ya que los datos cuentan con un proceso de curaduría previo en el que identificamos las especies que solo se encuentran en territorio nacional y aquellas dentro de alguna categoría de amenaza de acuerdo a la UICN.
# Extraemos la matriz de presencia-ausencia como data frame
paMa_df <- as.data.frame(paMa$Presence_and_Absence_Matrix)
#Filtramos las especies endémicas y amenazadas desde nuestra lista base
endemics <- filter(species, Endemic == "Yes")
threatened <- filter(species, IUCN.2023.1 %in% c("VU", "EN", "CR"))
# Validamos cuáles de las especies endémicas y amenazadas tienen polígono disponible
endemics_poly <- which(colnames(paMa_df) %in% endemics$Polygon.name)
threatened_poly <- which(colnames(paMa_df) %in% threatened$Polygon.name)
# Filtramos la matriz de presencia-ausencia a las especies endémicas y amenazadas
# --- Seleccionamos las columnas 1 y 2 que son lon y lat + los índices de las especies
paMa_df_endemics <- paMa_df[, c(1:2, endemics_poly)]
paMa_df_threatened <- paMa_df[, c(1:2, threatened_poly)]
# Sumamos el total de presencias para cada celda
# ---Como la matriz solo tiene 0/1, el total es el número de especies por píxel
paMa_df_endemics$Total <- rowSums(paMa_df_endemics[, 3:ncol(paMa_df_endemics)])
paMa_df_threatened$Total <- rowSums(paMa_df_threatened[, 3:ncol(paMa_df_threatened)])
# Renombramos las columnas de longitud y latitud
colnames(paMa_df_endemics)[1:2] <- c("x", "y")
colnames(paMa_df_threatened)[1:2] <- c("x", "y")
# Generamos una versión simplificada de la matriz solo con lon/lat y número de especies
paMa_df_endemics_s <- cbind.data.frame(x = paMa_df_endemics$x,
y = paMa_df_endemics$y,
Total = paMa_df_endemics$Total)
head(paMa_df_endemics_s)## x y Total
## 1 -83.75 13.75 0
## 2 -83.25 13.75 0
## 3 -83.75 13.25 0
## 4 -81.25 13.25 2
## 5 -83.75 12.75 0
## 6 -83.25 12.75 0paMa_df_threatened_s <- cbind.data.frame(x = paMa_df_threatened$x,
y = paMa_df_threatened$y,
Total = paMa_df_threatened$Total)
head(paMa_df_threatened_s)## x y Total
## 1 -83.75 13.75 0
## 2 -83.25 13.75 0
## 3 -83.75 13.25 0
## 4 -81.25 13.25 0
## 5 -83.75 12.75 0
## 6 -83.25 12.75 0# Convertimos los data frame de endémicas y amenazadas a raster
endemics_ras <- rasterFromXYZ(paMa_df_endemics_s[, 1:3])
threatened_ras <- rasterFromXYZ(paMa_df_threatened_s[, 1:3])
# Convertimos los raster de endémicas y amenazadas a un SpatRast autosuficiente
endemics_ras <- as(endemics_ras, "Raster") %>% as("SpatRaster")
endemics_ras <- terra::wrap(endemics_ras)
threatened_ras <- as(threatened_ras, "Raster") %>% as("SpatRaster")
threatened_ras <- terra::wrap(threatened_ras)
# Graficamos
plot(terra::unwrap(endemics_ras), col = pal)
plot(terra::unwrap(threatened_ras), col = pal)
5. Endemismo filogenético
Otra de las métricas interesantes que podemos calcular con los datos disponibles es el endemismo filogenético. Aquí comenzamos a integrar el componente espacial con el componente evolutivo. La idea de endemismo filogenético no es nueva (Rosauer et al. 2009) y combina dos métricas ampliamente utilizadas en la ecología filogenética de comunidades: la diversidad filogenética y el endemismo ponderado. No entraré en detalles, pero sí es importante aclarar que, el endemismo filogenético depende de dos factores principales: las longitudes de las ramas y el rango de distribución de las especies. De este modo, las comunidades que concentran especies con ramas más largas y rangos de distribución más pequeños tienden a presentar valores más altos.
Para calcular esta métrica solo necesitamos la matriz de presencia-ausencia y un árbol. Los árboles que se utilizan comúnmente para este tipo de análisis no necesariamente deben estar calibrados, aunque suele ser lo más habitual por razones de interpretación. Eso sí, deben estar enraizados y ser bifurcados. Sin embargo, como se mencionó antes, no contamos con un árbol completamente dicotómico que incluya todas nuestras especies de interés, por lo que estamos utilizando una distribución potencial de árboles. En este caso, se trata de árboles que fueron dicotomizados al azar para aquellos grupos sin información genética (Tonini et al. 2016). Dado que no tenemos uno sino muchos árboles con incertidumbre, esto implica un paso adicional: realizar el cálculo para cada uno y obtener un valor promedio.
# Creamos una copia de la matriz de presencia-ausencia (sin lon/lat)
paMa_df_s <- paMa_df[, -c(1:2)]
# Identificamos las especies con polígono que están en filogenia
# ---La lógica se basa en buscar especies que tengan datos en ambas columnas
commons <- species %>% filter(!is.na(Polygon.name) & !is.na(Phylo.name))
# Filtramos la matriz de presencia-ausencia a las especies comunes
# ---Usamos "Polygon.name" porque la matriz se construyó con base en los polígonos
paMa_df_s <- paMa_df_s %>% dplyr::select(any_of(commons$Polygon.name))
# Extraemos los nombres para cada polígono según como están en el árbol
# ---La idea es buscar los nombres de "Phylo.name" con base en los de "Polygon.name"
new_names <- species$Phylo.name[match(colnames(paMa_df_s), species$Polygon.name)]
# Reemplazamos los nombres de las columnas de la matriz de presencia-ausencia por los del árbol
colnames(paMa_df_s) <- new_names
# Podamos de nuevo los árboles pero solo a las especies con polígono
treepoly <- keep.tip(trees_pruned, colnames(paMa_df_s))
# Creamos una matriz vacía para almacenar los resultados
phyloen <- matrix(ncol = length(treepoly), nrow = nrow(paMa_df_s))
# Generamos un bucle para recorrer cada árbol
for(i in 1:length(treepoly)){
# Reordenamos las columnas de la matriz de presencia-ausencia
# ---Aunque es la misma matriz, el paquete que usaremos requiere una coincidencia
# ---perfecta en el órden de las especies entre la matriz y el árbol, y al ser
# ---árboles aleatorizados, cada uno puede tener las especies en una posición diferente
indices <- match(treepoly[[i]]$tip.label, colnames(paMa_df_s))
tmp.paMa <- paMa_df_s[, indices]
# Convertimos nuestra matriz a formato disperso (así lo pide el paquete)
tmp.paMa <- as(tmp.paMa, "sparseMatrix")
# Asignamos a cada pixel un número (sin esto, el paquete no funciona)
tmp.paMa@Dimnames[[1]] <- seq(1, nrow(tmp.paMa), by = 1)
# Calculamos el endemismo filogenético usando cada árbol y lo guardamos en la lista
phyloen[, i] <- phylo_endemism(tmp.paMa, treepoly[[i]])
}
# Calculamos el promedio de todos los valores de endemismo filogenético por pixel
# y asignamos de nuevo las coordenadas de cada celda desde la matriz de presencia-ausencia
phyloen_mean <- data.frame(paMa_df[, c(1:2)], phyloen = rowMeans(phyloen))
# Convertimos el data frame de endemismo filogenético a raster
phyloen_ras <- rasterFromXYZ(phyloen_mean[, 1:3])
# Convertimos el raster de endemismo filogenético a un SpatRast autosuficiente
phyloen_ras <- as(phyloen_ras, "Raster") %>% as("SpatRaster")
phyloen_ras <- terra::wrap(phyloen_ras)
# Graficamos
plot(terra::unwrap(phyloen_ras), col = pal)
6 Similitud filogenética
Para nuestro siguiente acto calcularemos la distancia filogenética promedio. Si bien esta no es una medida de diversidad filogenética en sentido estricto, la utilizaremos para hacer inferencias sobre qué tan parecidas o diferentes son nuestras comunidades desde un punto de vista evolutivo. Esto se debe a que la métrica más utilizada para medir diversidad filogenética, el índice de Faith (Faith 1992), es un surrogado filogenético de la diversidad α, por lo que, en términos espaciales, el patrón resultante sería idéntico.
Aquí, nuevamente, debemos considerar que estamos trabajando con un
conjunto de árboles, por lo que es necesario calcular las distancias
para cada uno y luego obtener un promedio. Existen muchos paquetes que
permiten calcular este tipo de métricas filogenéticas, pero
personalmente prefiero PhyloMeasures. Si bien es un paquete
antiguo y que ya no se actualiza, fue escrito en C++, lo que lo hace muy
eficiente y capaz de realizar cálculos para miles de comunidades en
cuestión de segundos.
Como en el cálculo del endemismo filogenético ya habíamos creado un bucle y filtrado las comunidades en función del árbol, aquí podemos reutilizar esos objetos.
# Creamos una matriz vacía para almacenar los resultados
mpd <- matrix(ncol = length(treepoly), nrow = nrow(paMa_df_s))
# Generamos un bucle para recorrer cada árbol
for(i in 1:length(treepoly)){
# Reordenamos las columnas de la matriz de presencia-ausencia
indices <- match(treepoly[[i]]$tip.label, colnames(paMa_df_s))
tmp.paMa <- paMa_df_s[, indices]
# Calculamos el MPD usando cada árbol y lo guardamos en la lista
mpd[, i] <- mpd.query(treepoly[[i]], tmp.paMa, reps = 999, standardize = FALSE)
}
# Calculamos el promedio de todos los valores de MPD por pixel y asignamos coordenadas
mpd_mean <- data.frame(paMa_df[, c(1:2)], mpd = rowMeans(mpd))
# Convertimos el data frame de MPD a raster
mpd_ras <- rasterFromXYZ(mpd_mean[, 1:3])
# Convertimos el raster de MPD a un SpatRast autosuficiente
mpd_ras <- as(mpd_ras, "Raster") %>% as("SpatRaster")
mpd_ras <- terra::wrap(mpd_ras)
# Graficamos
plot(terra::unwrap(mpd_ras), col = pal)
Nótese que, para el cálculo de la distancia filogenética, utilizamos
el argumento standardize = FALSE, lo que significa que no
estamos normalizando los valores. Normalizar, en este contexto, implica
transformar los datos a una distribución z con media cero y desviación
estándar uno. Esto es útil cuando queremos evaluar si el valor de una
comunidad en particular difiere de lo esperado por azar, ya que en una
distribución z, para un α = 0.05, los valores mayores a ±1.96, ya sean
positivos o negativos, se consideran significativos.
En otras palabras, la normalización facilita la interpretación de los resultados. En este caso, como estamos trabajando con miles de comunidades, no nos interesa evaluar la significancia estadística de cada una. Más adelante veremos un ejemplo en el que esto sí será relevante. Por otro lado, al no estandarizar los valores, los resultados se expresan en las mismas unidades que las longitudes de rama del árbol, en este caso, millones de años.
A simple vista, podemos notar que los valores del lado izquierdo del mapa parecen ser mayores que los del lado derecho. Este es el tipo de patrón que nos interesa explorar, ya que, a nivel biogeográfico, el lado izquierdo corresponde a grandes rasgos a la región transandina (donde están los Andes), mientras que el lado derecho corresponde a la región cisandina (donde está el Amazonas). La teoría nos dice que, la historia de las areas es también la historia de las especies, por lo que es posible que existan procesos eco-evolutivos detrás de este patrón. En este caso, realizaremos una prueba muy simple para evaluar si la mediana de ambas regiones difiere.
# Cargamos los polígonos de las regiones naturales de Colombia
regions <- st_read(paste0(getwd(), "/geo"), "Regiones_IGAC_Colombia", quiet = TRUE)
# Fusionamos los features según si pertenecen a la región cis- o trans-
regions_s <- regions %>% group_by(Portion) %>% summarise() %>% st_remove_holes()
# Enmascaramos el raster de MPD para extraer los valores de cada región
# ---Si el centroide del pixel no está dentro del polígono se elimina, por eso
# ---usamos el argumento touches = true, para prevenir este comportamiento
mpd_trans <- mask(terra::unwrap(mpd_ras), regions_s[2, ], touches = TRUE)
mpd_cis <- mask(terra::unwrap(mpd_ras), regions_s[1, ], touches = TRUE)
# Graficamos
plot(mpd_trans, col = pal)
plot(mpd_cis, col = pal)
# Extraemos los valores de los rasters
mpd_vals <- list(trans = values(mpd_trans, mat = FALSE, na.rm = TRUE),
cis = values(mpd_cis, mat = FALSE, na.rm = TRUE))
# Validamos la normalidad de los datos
# ---Si el valor p < 0.05 se rechaza la normalidad
# ---Cuando los datos son normales se usa ANOVA, cuando no, Kruskal-Wallis (no paramétrico)
lapply(mpd_vals, shapiro.test)## $trans
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.9649, p-value = 3.772e-05
##
##
## $cis
##
## Shapiro-Wilk normality test
##
## data: X[[i]]
## W = 0.96074, p-value = 1.416e-06# Calculamos un Kruskal-Wallis para comparar las medianas de los MPD de ambas regiones
kruskal.test(mpd_vals)##
## Kruskal-Wallis rank sum test
##
## data: mpd_vals
## Kruskal-Wallis chi-squared = 155.34, df = 1, p-value < 2.2e-16# Visualizamos las diferencias mediante un boxplot
boxplot(list(trans = values(mpd_trans, mat = FALSE, na.rm = TRUE),
cis = values(mpd_cis, mat = FALSE, na.rm = TRUE)))
Otra aclaración importante es que, la distancia que calculamos
anteriormente corresponde a la distancia filogenética promedio (MPD), la
cual se obtiene a partir de todos los pares posibles de especies dentro
de una comunidad. En cambio, existe otra métrica denominada distancia
media al taxón más cercano (MNTD), que solo considera un par por
especie, correspondiente al taxón evolutivamente más próximo. Para
calcular el MNTD en lugar del MPD, basta con reemplazar
mpd.query por mntd.query. Ambas métricas se
utilizan en contextos distintos y producen valores diferentes, por lo
que recomiendo revisar material externo para una explicación más
detallada, por ejemplo este libro: Phylogenetic
Ecology (Swenson 2019).
7 Mapas bivariados
Cuando tenemos datos de dos variables y queremos representarlas simultáneamente en el espacio geográfico, la mejor opción es utilizar mapas bivariados. En este tipo de representaciones se construye una matriz de colores bidimensional a partir de los rangos de ambas variables, por ejemplo usando cuantiles, que luego se utiliza para colorear áreas específicas. En este caso, emplearemos un mapa bivariado para mostrar de forma gráfica qué subregiones presentan un mayor número de especies y/o endemismos.
Como tenemos pocas regiones, recrearemos la tabla de datos directamente en R por simplicidad. Sin embargo, la extracción de esta información requiere algunos pasos adicionales. Para esto, utilicé de nuevo QGIS para calcular el número de polígonos que intersectan cada región en el caso de la diversidad total, así como el número de polígonos que se encuentran exclusivamente dentro de cada una para los endemismos. Si se preguntan por qué no hacerlo directo en R, es porque tengo una política, no hacer complicado lo que se puede hacer simple.
# Recreamos la tabla de datos para la diversidad y endemismo por subregión
bivariate <- data.frame(Subregion = c("Amazonas", "Andes", "Caribe", "Pacifico", "Orinoquia"),
Endemic = c(11, 128, 50, 35, 17),
Total = c(159, 357, 164, 171, 151))
# Unimos los datos por regiones con su geometría
bivariate <- regions %>% left_join(bivariate, by = "Subregion")
# Creamos las clases bivariadas
classes <- bi_class(bivariate,
x = Total,
y = Endemic,
style = "quantile",
dim = 3) # Número de particiones (en este caso, 3 cuantiles)
# Creamos una columna de colores para las etiquetas del mapa
# ---Esto es meramente estético, como Andes y Caribe quedan sobre fondo oscuro, se asigna un color claro
classes <- classes %>% mutate(color = ifelse(Subregion %in% c("Caribe", "Andes"), "white", "black"))
# Generamos el mapa bivariado (mucho de este código es solo carpientería para que se vea bien)
biplot <- ggplot(data = classes) +
geom_sf(aes(fill = bi_class), color = "black", size = 0.5, show.legend = FALSE) +
bi_scale_fill(pal = "DkViolet", dim = 3) +
geom_sf_text(
data = subset(classes, Subregion != "Pacifico"),
aes(label = Subregion, color = color),
size = 3) +
geom_sf_text(
data = subset(classes, Subregion == "Pacifico"),
aes(label = Subregion, color = color),
size = 3,
nudge_x = -1.8) +
scale_color_identity() +
theme_minimal() +
theme(
panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank(),
legend.position = "none")
# Generamos la leyenda del gráfico
legend <- bi_legend(pal = "DkViolet",
dim = 3,
xlab = "Riqueza",
ylab = "Endemismo",
size = 8)
# Unimos el mapa y la leyenda en un solo gráfico
ggdraw() +
draw_plot(biplot, 0, 0, 1, 1) +
draw_plot(legend, 0.05, 0.20, 0.25, 0.25)
Como en este análisis nuestras comunidades ya no son miles de píxeles, sino un número más manejable de regiones, podemos calcular la distancia filogenética promedio entre las especies de cada región en su versión estandarizada. Recordemos que la estandarización facilita la interpretación, ya que nos permite evaluar si el patrón observado difiere de lo esperado por azar. En este contexto, esto nos permite generar hipótesis sobre posibles reglas de ensamblaje de las comunidades, siempre que tengamos claro qué patrones esperar bajo distintos escenarios (Ochoa-Ochoa et al. 2020; Vásquez-Restrepo et al. 2023). Por ejemplo, una mayor dispersión de los organismos tiende a aumentar la dispersión filogenética, dando lugar a comunidades más disímiles, mientras que la diversificación in situ suele generar agrupamiento filogenético, es decir, comunidades más similares.
# Validamos que nuestro mapa de subregiones y polígonos tengan el mismo sistema de referencia
# ---Esta comparación debe devolver TRUE, de lo contrario, debemos estandarizar antes de continuar
st_crs(regions) == st_crs(IUCN)## [1] TRUE# Validamos y corregimos las geometrías por si hay alguna inválida
IUCN_sf <- st_as_sf(IUCN) %>% st_repair_geometry()
regions <- regions %>% st_repair_geometry()
# Intersectamos ambos objetos espaciales
region_species <- st_join(regions, IUCN_sf, left = FALSE)
# Generamos una matriz de presencia ausencia
paMa_regions <- region_species %>% st_drop_geometry() %>% dplyr::select(Subregion, BINOMIAL) %>% table()
paMa_regions <- as.data.frame(unclass(paMa_regions))
# Filtramos la matriz de presencia-ausencia a las especies comunes entre los polígonos y el árbol
# ---Usamos "Polygon.name" porque la matriz se construyó con base en los polígonos
paMa_regions <- paMa_regions %>% dplyr::select(any_of(commons$Polygon.name))
# Extraemos los nombres para cada polígono según como están en el árbol
# ---La idea es buscar los nombres de "Phylo.name" con base en los de "Polygon.name"
new_names <- species$Phylo.name[match(colnames(paMa_regions), species$Polygon.name)]
# Reemplazamos los nombres de las columnas de la matriz de presencia-ausencia por los del árbol
colnames(paMa_regions) <- new_names
# Podamos los árboles a las especies de esta nueva matriz de presencia-ausencia
# ---Como los polígonos para subregiones son continentales, vamos a perder las especies insulares
trees_subregion <- keep.tip(treepoly, new_names)
# Creamos una matriz vacía para almacenar los resultados
mpd_regions <- matrix(ncol = length(trees_subregion), nrow = nrow(paMa_regions))
# Generamos un bucle para recorrer cada árbol
for(i in 1:length(trees_subregion)){
# Reordenamos las columnas de la matriz de presencia-ausencia
indices <- match(trees_subregion[[i]]$tip.label, colnames(paMa_regions))
tmp.paMa <- paMa_regions[, indices]
# Calculamos el MPD usando cada árbol y lo guardamos en la lista
# ---Como tenemos pocas regiones, estandarizamos para comparar entre ellas
mpd_regions[, i] <- mpd.query(trees_subregion[[i]], tmp.paMa, reps = 999, standardize = TRUE)
}
# Calculamos el promedio de todos los valores de MPD por region y asignamos nombres
mpd_regions_mean <- data.frame(Region = rownames(paMa_regions), mpd = rowMeans(mpd_regions))
# Imprimimos el resultado
print(mpd_regions_mean)## Region mpd
## 1 Amazonas -2.4651741
## 2 Andes -0.8223060
## 3 Caribe 0.8188716
## 4 Orinoquia -4.1969976
## 5 Pacifico -0.4613242Con esto queda claro que solo el Amazonas y la Orinoquía presentan composiciones filogenéticas que difieren de lo esperado por azar, ya que sus valores superan el umbral de ±1.96. En particular, los valores negativos indican la presencia de agrupamiento filogenético.
8 Repensando las comunidades
Ahora, ¿qué pasaría si repensamos la idea tradicional de comunidad? Es decir, ¿qué tal si, en lugar de considerarlas como grupos de especies definidos por un contexto espacial, las entendemos como conjuntos de especies dentro de otro tipo de contexto? Por ejemplo, las categorías de la UICN. Hacer esta abstracción nos permite llevar los análisis un paso más allá, ya que podemos evaluar, por ejemplo, la contribución de las especies en cada categoría de amenaza a la historia evolutiva de nuestro ensamble global de especies.
Para este análisis calcularemos la diversidad filogenética dentro de cada categoría de la UICN utilizando dos métricas que mencionamos previamente: el índice de Faith y el MNTD. En el caso del índice de Faith, este no es más que la suma de las longitudes de rama que conectan a todas las especies de una comunidad hasta la raíz del árbol. Ahora bien, ¿por qué utilizar MNTD y no MPD en este contexto? Una diferencia clave es que el MNTD es más sensible a los patrones hacia la punta de las ramas, por lo que resulta más adecuado cuando analizamos estados de amenaza, que pueden interpretarse como atributos en el tiempo actual.
Diversidad filogenética promedio
# Filtramos las especies que están en filogenia
species_phylo <- species %>% filter(!is.na(Phylo.name))
# Generamos una pseudo matriz de presencia-ausencia para las categorías de la UICN
# ---La lógica es la misma, las comunidades en las filas y las especies en las columnas
paMa_IUCN <- table(species_phylo$IUCN.2023.1, species_phylo$Phylo.name)
paMa_IUCN <- unclass(paMa_IUCN)
# Creamos una matriz vacía para almacenar los resultados del índice de Faith
faith <- matrix(ncol = length(treepoly), nrow = nrow(paMa_IUCN))
# Generamos un bucle para calcular el índice de Faith usando cada árbol
for(i in 1:length(trees_pruned)){
# Calculamos la diversidad filogenética
tmp.PD <- tmp.PD <- pd(paMa_IUCN, trees_pruned[[i]])
# Guardamos la diversidad filogenética relativa (dividiendo sobre el número de especies)
# ---Esto lo hacemos porque Faith aumenta con la riqueza y las categorías son asimétricas
faith[, i] <- tmp.PD$PD / tmp.PD$SR
}
# Calculamos el promedio de todos los valores de Faith para los 100 árboles
faith_mean <- data.frame(Category = rownames(paMa_IUCN),
Mean = rowMeans(faith), # Promedio
SD = rowSds(faith)) # Desviación estandar
# Graficamos la diversidad filogenética promedio (mucho de este código es carpientería para que se vea bien)
faith_plot <- ggplot(data = faith_mean,
aes(x = factor(Category, levels = c("NE", "DD", "LC", "NT", "VU", "EN", "CR")),
y = Mean,
fill = Category)) +
geom_bar(stat = "identity") +
geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = 0.2, linewidth = 0.5) +
scale_fill_manual(values = c("#c52412", "#c3c3c3", "#f28533", "#51bc1d", "#e5e5e5", "#97c115", "#ffc90e")) +
theme_bw() +
theme(axis.text.x = element_text(size = 14, color = "black", margin = margin(t = 10)),
axis.text.y = element_text(size = 14, color = "black", margin = margin(r = 10)),
axis.ticks.length = unit(0.5, "cm"),
legend.position = "none",
strip.text.x = element_text(size = 14, margin = margin(0.9, 0, 0.9, 0, "cm")),
panel.background = element_blank(),
panel.grid = element_blank(),
panel.spacing = unit(1.5, "lines")) +
labs(x = NULL, y = NULL, title = NULL) +
scale_y_continuous(limits = c(0, 120), expand = c(0, 0), breaks = seq(0, 120, by = 20)) +
facet_wrap(~ "Diversidad filogenética promedio (Faith PD)")
print(faith_plot)
En la gráfica anterior podemos observar la contribución relativa promedio de las especies a cada categoría de la UICN, en términos de tiempo evolutivo. Es notable que las especies incluidas en alguna de las tres categorías de amenaza (VU, EN y CR) presentan una contribución alta. No obstante, este patrón podría cambiar en el futuro, dependiendo de cómo se clasifiquen las especies que aún no han sido evaluadas.
Distribución del estado de amenaza en la filogenia
# Creamos una matriz vacía para almacenar los resultados del MNTD
mntd <- matrix(ncol = length(treepoly), nrow = nrow(paMa_IUCN))
# Generamos un bucle para calcular el MNTD usando cada árbol
for(i in 1:length(trees_pruned)){
# Calculamos el mpd usando cada árbol y lo guardamos en la lista
# ---Como aquí tenemos menos comunidades, tiene sentido comparar, así que estandarizamos
mntd[, i] <- mntd.query(trees_pruned[[i]], paMa_IUCN, reps = 999, standardize = TRUE)
}
# Calculamos el promedio de todos los valores de MNTD para los 100 árboles
mntd_mean <- data.frame(Category = rownames(paMa_IUCN),
Mean = rowMeans(mntd), # Promedio
SD = rowSds(mntd)) # Desviación estandar
# Graficamos la distancia media al taxón más cercano
mntd_plot <- ggplot(data = mntd_mean,
aes(x = Mean,
y = factor(Category, levels = c("NE", "DD", "LC", "NT", "VU", "EN", "CR")),
fill = Category)) +
geom_errorbar(aes(xmin = Mean - SD, xmax = Mean + SD), width = 0.2, linewidth = 0.5) +
geom_point(size = 5, shape = 21, color = "black") +
geom_vline(xintercept = -1.96, size = 0.7, color = "black", linetype = "dashed") +
geom_vline(xintercept = 1.96, size = 0.7, color = "black", linetype = "dashed") +
scale_fill_manual(values = c("#c52412", "#c3c3c3", "#f28533", "#51bc1d", "#e5e5e5", "#97c115", "#ffc90e")) +
theme_bw() +
theme(axis.text.x = element_text(size = 14, color = "black", margin = margin(t = 10)),
axis.text.y = element_text(size = 14, color = "black", margin = margin(r = 10)),
axis.ticks.length = unit(0.5, "cm"),
legend.position = "none",
strip.text.x = element_text(size = 14, margin = margin(0.9, 0, 0.9, 0, "cm")),
panel.background = element_blank(),
panel.grid = element_blank(),
panel.spacing = unit(1.5, "lines")) +
labs(x = NULL, y = NULL, title = NULL) +
scale_x_continuous(limits = c(-3, 3), breaks = seq(-3, 3, 1), labels = seq(-3, 3, 1)) +
facet_wrap(~ "Distancia media al taxón más cercano (MNTD)")
print(mntd_plot)
Por otro lado, también es posible concluir que solo la categoría DD muestra un patrón de agrupamiento filogenético (valores negativos de MNTD) que difiere de lo esperado por azar, mientras que las demás categorías se ubican dentro de la región de no significancia. En términos simples, las categorías de amenaza se distribuyen de manera aleatoria a lo largo de la filogenia, con la excepción de las especies con datos deficientes, que tienden a agruparse en ciertos linajes. Este patrón también puede observarse fácilmente en el árbol filogenético.
# Elegimos un árbol al azar para representar la figura
random_tree <- trees_pruned[[26]]
# Convertimos el árbol a formato tabla (necesario para poder mapear los nodos)
random_tree_df <- as_tibble(random_tree)
# Añadimos la categoría de la UICN a cada especie en el árbol-tabla
random_tree_df <- random_tree_df %>% left_join(species_phylo[, c(10, 20)], by = c("label" = "Phylo.name"))
# Creamos un vector de colores según las categorías de la UICN
colores_IUCN <- c("NE" = "#e5e5e5",
"DD" = "#c3c3c3",
"LC" = "#51bc1d",
"NT" = "#97c115",
"VU" = "#ffc90e",
"EN" = "#f28533",
"CR" = "#c52412")
# Añadimos los colores por categoría al árbol-tabla
random_tree_df <- random_tree_df %>% mutate(color = colores_IUCN[IUCN.2023.1])
# Asignamos un color gris a las ramas internas del árbol
random_tree_df$color[which(is.na(random_tree_df$color))] <- "#bebebe"
# Graficamos el árbol y asignamos los colores por rama
IUCN_tree <- ggtree(random_tree, layout = "fan", open.angle = 10, size = 0.5) %<+% random_tree_df +
aes(color = I(color))
print (IUCN_tree)
Aquí lo tenemos, una forma visual de observar cómo se distribuyen los estados de amenaza a lo largo de una filogenia. Sin embargo, esta gráfica es bastante básica en cuanto a la información que presenta, ya que, debido al número de especies, incluir los nombres de cada una la haría ilegible. Por esta razón, resulta más adecuado, desde un punto de vista estilístico, utilizar etiquetas de grupos mayores, como familias.
Aunque este proceso también puede realizarse en R, hacerlo es programáticamente más complejo. Una alternativa más práctica es utilizar un software externo como TreeViewer. Este programa ofrece opciones avanzadas de edición mediante una interfaz gráfica bastante intuitiva. En este caso, solo sería necesario exportar la tabla que acabamos de crear y añadir columnas con el o los grupos que queremos resaltar, algo relativamente fácil de hacer a partir de los nombres de las especies.
Finalmente, podemos respaldar nuestra conclusión sobre la dispersión de las categorías de amenaza en la filogenia utilizando una medida de señal filogenética. En este caso, dado que trabajamos con caracteres discretos y cualitativos, métricas tradicionales para datos continuos como el lambda de Pagel o el K de Blomberg no son adecuadas. Por ello, utilizaremos el estadístico D de Fritz-Purvis (Fritz & Purvis 2010).
# Creamos una variable binaria para amenazada (1) y no amenazada (0)
# ---Esto es necesario porque el estadístico D solo analiza caracteres binarios
# ---La lógica del ifelse es simple, si se cumple la condición asigna 1, si no, 0
species_phylo$IUCN.bin <- ifelse(species_phylo$IUCN.2023.1 %in% c("VU", "EN", "CR"), 1, 0)
# Creamos una lista vacía para almacenar los resultados
D <- list()
# Generamos un bucle para calcular el estadístico D usando cada árbol
for(i in 1:length(trees_pruned)){
D[[i]] <- phylo.d(data = species_phylo, phy = trees_pruned[[i]], names.col = Phylo.name, binvar = IUCN.bin)
}
# Calculamos el promedio del estadístico D a partir de todas las iteraciones
mean(sapply(D, function(x) x$DEstimate), na.rm = TRUE)## [1] 0.7897616# Extraemos el rango de los valores para el estadístico D
d_vals <- sapply(D, function(x) x$DEstimate)
# Visualizamos el rango de los valores para el estadístico D
hist(d_vals, main = NULL, xlab = "Purvis D")
Este estadístico tiene una interpretación directa y sencilla: valores cercanos a cero indican la presencia de señal filogenética, es decir, una evolución no aleatoria, mientras que valores cercanos a uno indican ausencia de señal, consistente con una evolución aleatoria. En este caso, como podemos observar, el valor promedio se aleja de cero, con una tendencia a valores altos (> 0.7).
9 Priorizando especies desde el punto de vista evolutivo
Finalmente, quisiera cerrar este tutorial con un último análisis, uno que nos permite generar una lista de prioridades de conservación desde una perspectiva evolutiva. Tradicionalmente, las especies se priorizan con base en su categoría de amenaza, la cual se fundamenta en una serie de criterios geográficos y poblacionales (usualmente más geográficos que poblacionales). Estos criterios nos hablan de cierta “rareza” de las especies y de su vulnerabilidad a la extinción. Sin embargo, pocas veces nos preguntamos cuál es su contribución a la historia evolutiva.
A partir de los análisis anteriores ya hemos visto qué áreas o categorías contribuyen más o menos a nivel filogenético, pero ahora podemos ir un paso más allá y poner nombre a las especies que generan esas contribuciones. Para ello utilizaremos una metodología conocida como EDGE (Gumbs et al. 2023), un enfoque de valor y riesgo que combina la categoría de amenaza actual de las especies con su historia evolutiva. Una vez más, aquí me limitaré a mostrar el procedimiento para calcular los valores EDGE, la teoría detrás del método y su interpretación quedan como tarea para el lector.
# Cargamos las funciones auxiliares para EDGE2
# ---No existe un paquete para calcular esta métrica, pero los autores disponibilizaron una
# ---serie de funciones que nos facilitan la vida, es como cargar una librería, pero local
source("auxfun/EDGE2.R")
# Extraemos los datos de especies y categoría UICN
spe2EDGE <- data.frame(Species = species_phylo$Phylo.name,
IUCN = species_phylo$IUCN.2023.1)
# Cargamos la distribución de probabilidades de extinción por categoría
# ---Para hacerlo más rápido, vamos a cargar una versión precompilada
# ---También se puede generar un nuevo conjunto de valores usando: GE2 <- GE.2.calc(pext.vals)
# ---Para entender la idea detrás de esto véase (Gumbs et al. 2023)
GE2 <- read.csv("auxfun/GE2.csv")
# Generamos para cada especie en cada árbol una serie de probabilidades de extinción
set.seed(955)
pext <- random_GE2(spe2EDGE, GE2, 100)
# Creamos una lista vacía para almacenar los resultados
EDGE2 <- list()
# Generamos un bucle para calcular EDGE usando cada árbol y sus respectivas probabilidades
for(i in 1:length(trees_pruned)){
EDGE2 <- append(EDGE2, EDGE2_mod(trees_pruned[[i]], pext[, c(1, (2 + i))])[1])
}
# Unimos todos los resultados en una sola tabla simplificada
EDGE2_s <- EDGE2 %>% map(dplyr::select, "Species", "EDGE") %>% reduce(left_join, by = "Species")
# Renombramos las columnas de la tabla simplificada
colnames(EDGE2_s)[2:ncol(EDGE2_s)] <- paste("EDGE", 1:(ncol(EDGE2_s) - 1), sep = ".")
# Calculamos las medianas de los valores EDGE
# ---Este paso es necesario para poder asignar categorías EDGE (véase Gumbs et al.)
thresholds <- colMedians(as.matrix(EDGE2_s[, -1]), na.rm = TRUE)
# Creamos un dataframe binario comparando la matriz completa contra los thresholds
# ---Acá, binarizar se traduce como EDGE (1) o no EDGE (0)
EDGE2_bi <- as.data.frame(sweep(as.matrix(EDGE2_s[, -1]), 2, thresholds, ">") * 1)
EDGE2_bi$Species <- EDGE2_s$Species # Reasignamos las especies
EDGE2_bi$SUM <- rowSums(EDGE2_bi[, -ncol(EDGE2_bi)]) # Sumamos filas
# Creamos una lista final y clasificamos
EDGE_list <- data.frame(Species = EDGE2_s$Species,
EDGE.median = rowMedians(as.matrix(EDGE2_s[, -1]), na.rm = TRUE),
times = EDGE2_bi$SUM) # Número de veces que la especie se clasificó EDGE
# Añadimos de nuevo las categorías de la UICN a nuestra lista final
EDGE_list <- left_join(EDGE_list, species_phylo[, c(20, 10)], by = c("Species" = "Phylo.name"))
# Renombramos la columna UICN
colnames(EDGE_list)[4] <- "IUCN"
# Clasificamos las especies según los grupos propuestos en Gumbs et al. (2023)
EDGE_list <- EDGE_list %>% mutate(Group = case_when(times < 80 ~ "Not EDGE",
times >= 80 & times < 95 & IUCN %in% c("VU", "EN", "CR") ~ "Borderline",
times >= 80 & times < 95 ~ "Not EDGE",
times >= 95 & IUCN %in% c("VU", "EN", "CR") ~ "EDGE",
times >= 95 & IUCN %in% c("NE", "DD") ~ "Research List",
times >= 95 & IUCN %in% c("LC", "NT") ~ "Watch List"))
# Renombramos las especies de la lista final con su nombre válido
EDGE_list$Species <- species_phylo$Species[match(EDGE_list$Species, species_phylo$Phylo.name)]
# Visualizamos las primeras 25 especies EDGE y/o Borderline
toFilter <- c("EDGE", "Borderline")
head(EDGE_list %>% filter(Group %in% toFilter) %>% arrange(match(Group, toFilter), desc(EDGE.median)), 25)## Species EDGE.median times IUCN Group
## 1 Sphaerodactylus scapularis 15.756773 100 EN EDGE
## 2 Lepidoblepharis miyatai 12.265058 98 CR EDGE
## 3 Emmochliophis miops 11.242819 100 CR EDGE
## 4 Coniophanes andresensis 9.193017 100 CR EDGE
## 5 Anolis calimae 7.690679 98 VU EDGE
## 6 Anolis ruizii 7.063515 97 EN EDGE
## 7 Alopoglossus lehmanni 6.671974 95 CR EDGE
## 8 Lepidoblepharis williamsi 6.604680 99 EN EDGE
## 9 Alopoglossus danieli 6.058519 99 CR EDGE
## 10 Micrurus medemi 5.447006 97 CR EDGE
## 11 Synophis plectovertebralis 4.709375 100 CR EDGE
## 12 Enyalioides groi 4.704739 100 EN EDGE
## 13 Enyalioides oshaughnessyi 3.902355 100 VU EDGE
## 14 Corallus blombergi 3.559714 96 EN EDGE
## 15 Bothrocophias campbelli 3.295269 99 VU EDGE
## 16 Synophis bicolor 3.231079 97 EN EDGE
## 17 Riama columbiana 3.158061 98 EN EDGE
## 18 Enyalioides annularis 2.991262 100 VU EDGE
## 19 Riama simotera 2.743679 97 EN EDGE
## 20 Dendrophidion boshelli 3.115995 91 CR Borderline
## 21 Dipsas ellipsifera 2.434029 93 EN Borderline
## 22 Anadia pamplonensis 1.770867 81 EN Borderline
## 23 Saphenophis sneiderni 1.763351 91 EN Borderline
## 24 Andinosaura laevis 1.479432 91 VU Borderline
## 25 Micrurus sangilensis 1.168866 81 VU Borderline