Regla de Rapoport y otras patrones

Axel Arango & Fabricio Villalobos

Regla de Rapoport

“Relación positiva entre latitud y tamaño de área de distribución” G. Stevens

En este ejercicio se probará la regla de Rapoport en un set de más de 6000 registros para especies de la familia Icteridae (Aves: Passeriformes) y se compararán los procedimientos utilizando dos paqueterías: letsR y EcoPhylo Mapper epm

epm utiliza los paquetes sf y terra para análisis espaciales, lo cuál será útil cuándo rgdal y raster se hayan descontinuado (a partir de ya! junio 2023)

Paquetes necesarios:

library(epm)
library(rgdal)
library(sf)
library(caper)
library(here)
library(geiger)
library(letsR)
library(ggplot2)
library(dplyr)
library(spdep)
library(spatialreg)

Cargar los datos de registros de Icteridae (obtenidos de GBIF)

ict<-read.csv(here("data","ict.csv"),header=T)

Limpiar los registros, asegurándonos que no haya ningún NA ni registros repetidos:

ict<-ict %>% 
  na.omit %>% 
  unique %>% 
  rename(taxon=sp,x=longitude,y=latitude)

Finalmente, cambiaremos los espacios por guiones bajos en los nombres de las especies. Esto con el fin de que puedan ser reconocidas al momento de utilizar alguna filogenía

ict<-ict %>% 
  mutate(taxon=gsub(" ","_",taxon))
head(ict)

                taxon        x       y
1     Icterus_wagleri -96.8249 17.1452
2        Cacicus_cela -80.7952 -1.6779
3 Icterus_graduacauda -96.8522 17.1685
4   Icterus_parisorum -96.8249 17.1452
5     Sturnella_magna -96.7706 17.1733
6     Icterus_wagleri -96.8522 17.1685

Habiendo cargado y limpiado los registros, podemos obtener los atributos de interés para evaluar la regla de Rapoport en esta familia de aves. Para esto necesitamos saber el área de distribución de las especies y su ubicación geográfica.

Usando letsR, primero, creamos una matriz de presencia-ausencia

presab<-lets.presab.points(cbind(ict$x,ict$y),ict$taxon)

Después, obtenemos el área de distribución de las especies

ranges<-lets.rangesize(presab)
head(ranges)

                     Range_size
Agelaioides_badius           36
Agelaius_phoeniceus         269
Agelaius_tricolor            12
Agelaius_xanthomus            1
Agelasticus_cyanopus          4
Agelasticus_thilius          15

Finalmente, necesitamos la ubicación de las especies.

Afortunadamente, letsR posee una práctica función que permite obtener el punto medio de la ubicación geográfica de las especies

midp<-lets.midpoint(presab)
head(midp)

               Species          x         y
1   Agelaioides_badius  -64.90127 -31.71035
2  Agelaius_phoeniceus  -97.27983  30.32207
3    Agelaius_tricolor -121.17026  36.84573
4   Agelaius_xanthomus  -67.50000  17.50000
5 Agelasticus_cyanopus  -57.50000 -28.50118
6  Agelasticus_thilius  -58.16536 -34.17206

Ahora que tenemos los datos, podemos poner a prueba la relación esperada por la Regla de Rapoport

lm(log(ranges)~abs(midp$y)) %>% 
summary


Call:
lm(formula = log(ranges) ~ abs(midp$y))

Residuals:
    Min      1Q  Median      3Q     Max 
-2.8064 -1.0526  0.2179  0.9897  2.7981 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.14109    0.29532   3.864 0.000235 ***
abs(midp$y)  0.05460    0.01252   4.361 4.07e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.233 on 75 degrees of freedom
Multiple R-squared:  0.2023,    Adjusted R-squared:  0.1916 
F-statistic: 19.01 on 1 and 75 DF,  p-value: 4.074e-05

La relación entre el área de distribución y la latitud es significativamente postiva en Icteridae ahora observemosla de manera más gráfica:

dfram<-data.frame(midp,ranges);names(dfram)<-c("sp","x","y","ranges")
p<-ggplot(dfram,aes(abs(y),log(ranges)))+
            geom_point(color="black")+
            xlab("Latitud absoluta")+
            ylab("Log (tamaño área)")+
            geom_smooth(method= lm , color="red", fill="#69b3a2", se=TRUE)+
            theme_classic()

plot(p)

Ahora utilicemos el paquete epm para tratar de hacer lo mismo. epm no puede utilizar sòlo los puntos de longitud y latitud de nuestros registros, primero hay que transformar esos puntos en un “spatial feature”

sp_ict<-st_as_sf(ict, coords = c("x","y"), crs= st_crs("+datum=WGS84 +proj=longlat"))
head(sp_ict)

Simple feature collection with 6 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -96.8522 ymin: -1.6779 xmax: -80.7952 ymax: 17.1733
Geodetic CRS:  +datum=WGS84 +proj=longlat
                taxon                 geometry
1     Icterus_wagleri POINT (-96.8249 17.1452)
2        Cacicus_cela POINT (-80.7952 -1.6779)
3 Icterus_graduacauda POINT (-96.8522 17.1685)
4   Icterus_parisorum POINT (-96.8249 17.1452)
5     Sturnella_magna POINT (-96.7706 17.1733)
6     Icterus_wagleri POINT (-96.8522 17.1685)

Ahora que nuestros registros son un objeto espacial, ya podemos utilizar epm para crear un grid con nuestras especies en el espacio:

grix<-createEPMgrid(sp_ict,resolution = 1)

Con este grid, ya podemos obtener el área de distribución de las especies, dónde cell es área de distribución

Cellcount<-data.frame(grix$cellCount)
names(Cellcount)<-"cell"
head(Cellcount)

                     cell
Agelaioides_badius     36
Agelaius_phoeniceus   287
Agelaius_tricolor      12
Agelaius_xanthomus      1
Agelasticus_cyanopus    4
Agelasticus_thilius    15

¿Qué tanto se parece la distribución calculada con letsR a la calculada con epm?

par(mfrow=c(1,2))
hist(dfram$ranges, breaks=15, main ="LetsR", xlab="área de distribución",ylab="especies")
hist(Cellcount$cell, breaks =15, main= "EPM", xlab ="área de distribución",ylab="especies")

Para la ubicación geográfica de las especies, en el caso de epm no se cuenta con una función especifica para obtener los puntos medios de la distribución de las especies.

Por lo tanto, tenemos que cargar una función personalizada:

load(here("data","st_midpoints.R"))

. . . Esta función utiliza la nueva paquetería sf para crear un polígono para cada especie con todas las coordenadas de nuestros registros, para posteriormente calcular el centroide de dicho polígono, obteniendo así los puntos medios de la ubicación geográfica de nuestras especies.

midst<-st_midpoints(ict)
head(midst)

                     sp          X         Y
55   Agelaioides_badius  -53.18129 -22.92205
18  Agelaius_phoeniceus  -97.27558  30.89165
44    Agelaius_tricolor -119.04709  35.69512
73   Agelaius_xanthomus  -67.11699  17.95263
15 Agelasticus_cyanopus  -53.72550 -21.23416
26  Agelasticus_thilius  -65.57207 -26.82030

Ahora ya tenemos los datos necesarios para evaluar la relación esperada por la Regla de Rapoport

lm(log(Cellcount$cell)~abs(midst$Y)) %>% 
  summary


Call:
lm(formula = log(Cellcount$cell) ~ abs(midst$Y))

Residuals:
     Min       1Q   Median       3Q      Max 
-2.98520 -0.72200  0.07052  0.86243  2.81757 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.19213    0.28009   4.256 5.94e-05 ***
abs(midst$Y)  0.05903    0.01278   4.620 1.56e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.212 on 75 degrees of freedom
Multiple R-squared:  0.2216,    Adjusted R-squared:  0.2112 
F-statistic: 21.35 on 1 and 75 DF,  p-value: 1.56e-05

El resultado es una relación significativamente postiva en Icteridae entre el área de distribución y la latitud.
Esta relación es muy parecida a la obtenida utilizando el procedimiento letsR. Las diferencias observadas pueden deberse a 1) el cálculo de los grids, pues letsR utiliza celdas y nosotros calculamos hexagonos con epm o 2) el cálculo de los puntos medios, pues para epm dicho cálculo es más o menos artesanal.
Ahora que hemos evaluado la regla de Rapoport utilizando dos métodos y paqueterías distintas, podemos pensar en el efecto del componente evolutivo sobre dicha relación.

Para esto, primero debemos cargar una filogenía para nuestras especies. En este caso utilizaremos una filogenía para Icteridae, cortada de la filogenía para Emberizoidea, usada en Arango et al., 2022, modificada y estandarizada de Barker et al., 2015.

tre<-read.tree(here("data","Icteridae_tree.txt"))
tre


Phylogenetic tree with 100 tips and 99 internal nodes.

Tip labels:
  Icterus_melanopsis, Icterus_northropi, Icterus_laudabilis, Icterus_dominicensis, Icterus_cayanensis, Icterus_pyrrhopterus, ...
Node labels:
  761, 762, 763, 764, 765, 766, ...

Rooted; includes branch lengths.

Ahora crearemos un data.frame el cual contenga el área de distribución y los puntos medios para cada especie. Para motivos del ejercicio, haremos uso de los resultados obtenidos con epm. ¡Pero ustedes pueden hacerlo con los resultados obtenidos con letsR y podemos comparar las relaciones finales!

rapp<-midst %>% 
  cbind(Cellcount) %>% 
  as.data.frame() %>% 
  rename(x=X,y=Y,area=cell)

head(rapp)

                     sp          x         y area
55   Agelaioides_badius  -53.18129 -22.92205   36
18  Agelaius_phoeniceus  -97.27558  30.89165  287
44    Agelaius_tricolor -119.04709  35.69512   12
73   Agelaius_xanthomus  -67.11699  17.95263    1
15 Agelasticus_cyanopus  -53.72550 -21.23416    4
26  Agelasticus_thilius  -65.57207 -26.82030   15

Una vez que tenemos nuestra tabla de datos, creamos un objeto comparative.data, el cual juntará nuestros datos geográficos con nuestros datos filogenéticos

compict<-comparative.data(tre,data =rapp,names.col = "sp",vcv = T)

Ahora ya podemos ajustar un modelo de regresión considerando las relaciones filogenéticas de nuestras especies utilizando un PGLS (Phylogenetic Generalized Least Squares):

pgls(log(area)~abs(y),data=compict,lambda="ML") %>% 
summary()


Call:
pgls(formula = log(area) ~ abs(y), data = compict, lambda = "ML")

Residuals:
     Min       1Q   Median       3Q      Max 
-1.07514 -0.24240  0.00397  0.16635  0.76134 

Branch length transformations:

kappa  [Fix]  : 1.000
lambda [ ML]  : 0.000
   lower bound : 0.000, p = 1    
   upper bound : 1.000, p = 4.5519e-14
   95.0% CI   : (NA, 0.651)
delta  [Fix]  : 1.000

Coefficients:
            Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 1.274120   0.296186  4.3018 5.551e-05 ***
abs(y)      0.057690   0.013493  4.2755 6.093e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3893 on 68 degrees of freedom
Multiple R-squared: 0.2119, Adjusted R-squared: 0.2003 
F-statistic: 18.28 on 1 and 68 DF,  p-value: 6.093e-05

¡La relación se mantiene!

Ahora, si bien es mucho más fácil obtener estos resultados utlizando letsR, epm hace uso de las nuevas paqueterías espaciales, además de servir para mapear otros componentes de diversidad y macroevolución.
Por ejemplo, podemos obtener la diversidad filogenética (PD), las tasas de diversificación (DR) y los endemismos filogenéticos de cada comunidad (o celda) y responder preguntas que involucren métricas de este tipo y la geografía.

Para esto, primero debemos agregarle una filogenía a nuestro grid creado con epm

grix<-epm::addPhylo(grix,tre)

Ya con la filogenía incorporada a nuestro grid de epm, calculamos para cada celda PD:

pd<-gridMetrics(grix,"pd")

Además, podemos calcular otras métricas ‘evolutivas’ como la tasa de diversificación y los endemismos filogenéticos

Tasa de especiación

dr<-gridMetrics(grix,"DR")

Endemismo filogenético

phyendm<-gridMetrics(grix,"phyloWeightedEndemism")

Veamos estas métricas mapeadas en el espacio geográfico

par(mfrow=c(1,3))
plot(pd,use_tmap=F, legend=F)
addLegend(pd,location="left",label = "PD")
plot(dr,use_tmap=F, legend=F)
addLegend(dr,location="left",label = "DR")
plot(phyendm,use_tmap=F, legend=F)
addLegend(phyendm,location="left",label = "Phyloendemisms")

Con estos datos, podemos hacernos preguntas cómo: ¿Las comunidades que poseen mayores tasas de diversificación tienen una mayor diversidad filogenética?

summary(lm(pd$grid$pd~dr$grid$DR))


Call:
lm(formula = pd$grid$pd ~ dr$grid$DR)

Residuals:
    Min      1Q  Median      3Q     Max 
-11.580  -8.297  -3.003   5.778  48.604 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  17.0500     0.9138  18.658   <2e-16 ***
dr$grid$DR    3.9509     1.8877   2.093   0.0367 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.79 on 738 degrees of freedom
  (27 observations deleted due to missingness)
Multiple R-squared:  0.005901,  Adjusted R-squared:  0.004554 
F-statistic: 4.381 on 1 and 738 DF,  p-value: 0.03669

plot(dr$grid$DR,pd$grid$pd,xlab="Tasa de diversificación",ylab="Diversidad filogenética", pch=16)

¿Los endemismos filognéticos siguen un patrón de gradiente latitudinal?

Obtengamos las coordenadas de cada celda y su valor de filoendemismos:

coordsx<-grix$grid$gridTemplate %>% 
  st_centroid %>% 
  st_coordinates %>% 
  as.data.frame
endmcoords<-phyendm$grid$phyloWeightedEndemism %>% 
  cbind(coordsx) %>% 
  as.data.frame %>% 
  na.omit 
names(endmcoords)<-c("phyloendemism","x","y")
head(endmcoords)

  phyloendemism         x        y
1    0.05185643 -132.1468 54.27687
2    0.05185643 -124.6468 43.01901
3    0.07705955 -124.1468 40.42101
4    0.07705955 -124.1468 43.88501
5    0.07242206 -124.1468 45.61700
6    0.07242206 -124.1468 49.08097

Ahora, pongamos a prueba la pregunta:

lm(log(phyloendemism)~abs(y),data= endmcoords) %>% 
  summary


Call:
lm(formula = log(phyloendemism) ~ abs(y), data = endmcoords)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.1211 -0.7451 -0.0409  0.6642  3.3798 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.515137   0.099321  -5.187 2.77e-07 ***
abs(y)      -0.049390   0.003193 -15.467  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.03 on 738 degrees of freedom
Multiple R-squared:  0.2448,    Adjusted R-squared:  0.2438 
F-statistic: 239.2 on 1 and 738 DF,  p-value: < 2.2e-16

Parece que los endemismos disminuyen al alejarse del ecuador, siguiendo un gradiente latitudinal. Observémoslo:

p2<-ggplot(endmcoords,aes(abs(y),log(phyloendemism)))+
            geom_point(color="black")+
            xlab("Latitud absoluta")+
            ylab("log (Filoendemismos)")+
            geom_smooth(method= lm , color="red", fill="#69b3a2", se=TRUE)+
            theme_classic()

plot(p2)

¿Qué podemos concluir de esto?